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X.  INTRODUCTION 


This  report  has  been  organized  into  six  main  sections. 

Section  1  introduces  and  discusses  the  purpose  of  this  report, 
the  research  objectives,  and  pertinent  background  material. 
Sections  2  through  4,  discuss  the  sources  of  computational  flow 
analysis  error,  suggest  a  procedure  (developed  during  the  course 
of  this  investigation)  for  monitoring  and  locating  these  errors, 
and  demonstrate  the  error  monitor  procedure  for  a  configuration 
of  engineering  interest. 

The  contract  investigation  raised  several  issues  regarding 
the  usage  and  generality  of  the  error  monitor.  These  issues  were 
addressed  in  a  research  effort  conducted  using  Boeing  Advanced 
Systems  (BAS)  Independent  Research  and  Development  (IRAD)  funds. 
The  results  of  these  research  studies  are  summarized  in  Section  5 
and  are  included  in  their  entirety  in  the  Appendices. 

The  last  section.  Section  6,  summarizes  the  status  of  the 
research.  Significant  accomplishments,  conclusions,  and 
recommendations  are  enumerated  in  this  section. 

The  reader  with  limited  time  may  wish  to  review  the 
Introduction  and  Research  Status  sections  first  (Sections  1  and 
6).  Sections  2  through  5  provide  a  more  detailed  discussion. 

The  Appendices  contain  the  detailed  technical  information  useful 
as  a  starting  point  for  further  investigations  of  the  error 
monitors  and  procedures  discussed  in  the  body  of  this  report. 


1.1  Purpose  of  Final  Report 

This  report  summarizes  the  results  of  work  performed  under 
the  contract  titled  "Error  Norm  Guided  Flow  Analysis  of  Shock- 
Wave/Boundary-Layer  Interactions"  contract  number  F49620-85-C- 
0126  and  a  related  follow-on  IRAD  study.  Previous  reports  by 
Forester  (1986,  2/87,  and  8/87)  have  documented  preliminary 
results  and  progress  of  the  contract  effort.  Pertinent 
information  from  these  documents  is  also  summarized  in  this,  the 
final  contract  report.  The  purpose  of  the  contract  was  to 
develop  a  grid  selection  and  smoothing  level  guide  for  accurately 
applying  Navier-Stokes  (NS)  analysis  to  flow  fields  of 
engineering  interest.  The  guide  indicates  grid  resolution 
problems  and  is  computed  from  data  readily  available  in  NS 
solutions.  The  grid  guide  was  formatted  for  graphical  output. 

An  inlet  design  with  a  shock  control  aperture  was  chosen  to 
demonstrate  the  concept. 

The  contract  results  were  encouraging,  but  raised  several 
questions  on  the  usage  and  generality  of  the  grid  selection 
guide.  These  questions  included  the  usefulness  of  the  guide  for 
algorithms  other  than  the  MacCormack  explicit  (used  for  the 
contract  study),  the  relationship  between  the  guide  and  solution 
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acc'-  acy,  and  the  usefulness  of  the  guide  for  a  wide  range  of 
|  grid  levels.  A  follow-on  IRAD  study  was  conducted  to  address 

these  questions.  In  the  follow-on  study,  additional  grid 
refinement  work  was  performed,  the  generality  of  the  error 
monitor  and  extensions  to  other  algorithms  were  investigated,  and 
the  relationship  between  the  error  monitor  magnitude  and 
truncation  error  levels  was  studied. 

I 

1.2  Research  Objectives 

The  overall  contract  objective  was  to  explore  grid  adjustment 
guides  for  improving  the  accuracy  of  flow  field  analyses.  To 
S  accomplish  this  the  following  specific  objectives  were 

1  identified: 

o  Develop  and  evaluate  error  norms  as  guides  to  grid 
selection  and  smoothing  level. 

o  Define  and  implement  these  error  norms  in  a  two-dimensional 
(2D)  NS  code. 

c  Demonstrate  and  evaluate  error  norm  usage  with  a  NS 
analysis  of  a  2D  inlet. 

This  work  was  accomplished  using  MacCormack's  explicit  predictor- 
I  !•  corrector  method  with  explicitly  added  numerical  smoothing. 


1.3  Background 

Computational  fluid  dynamic  ( CFD )  tools  are  needed  to  provide 
flow  field  properties  for  the  design  of  aircraft  components.  For 
three-dimensional  analyses,  Paynter  and  Chen,  and  Anderson  have 
made  considerable  progress.  Three  dimensional  CFD  analyses  are 
currently  used  in  the  selection  and  design  of  test  configurations 
and  full-scale  configurations.  Analysis  also  enhances  the  degree 
of  integration  of  the  propulsion  system  with  flight  controls, 
aircraft  aerodynamics,  and  structural  requirements.  The  proper 
integration  of  these  elements  is  the  key  to  high  performance 
aircraft  design. 

When  CFD  is  used  to  provide  flow  properties  for  a  flow  domain 
of  interest,  algebraic  approximations  to  the  NS  equations  (or  a 
subset)  are  solved  for  the  values  of  the  dependent  variables  at 
each  point  on  a  grid  of  points  distributed  through  the  domain. 
Available  computer  speed  and  storage  limit  the  number  of  points 
at  which  the  algebraic  approximations  to  the  NS  equations  (Finite 
difference  equations  or  FDEs)  can  be  solved  (for  a  given  analysis 
domain).  Details  of  the  flow  in  the  domain  are  unknown  prior  to 
obtaining  a  solution,  and  local  regions  with  high  gradients  in 
flow  properties  (such  as  shocks  or  shear  layers)  can  occur  and  be 
of  critical  importance  to  the  performance  of  an  aircraft 
component.  A  coarse  grid  in  a  region  of  high  local  flow  gradient 
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can  inadequately  resolve  these  gradients  and  result  in  a  high 
local  solution  error,  manifested  as  an  oscillation  in  flow 
properties  near  the  high  gradient  region  (dispersive  error). 
Dispersive  error  can  be  controlled  to  some  degree  at  least 
through  the  use  of  artificial  dissipation  (smoothing)  in  regions 
of  high  flow  gradient  and  coarse  mesh.  The  selection  of  a  grid 
point  distribution  and  smoothing  levels  to  achieve  a  desired 
level  of  solution  accuracy  are  recognized  as  major  problems  in 
the  use  of  CFD  for  aircraft  design  support. 

Flow  field  solution  accuracy  is  a  strong  function  of  grid 
density,  grid  placement  and  smoothing  levels.  Methods  are  needed 
for  controlling  solution  errors  due  to  grid  and  for  selecting 
numerical  smoothing  levels.  In  general,  simple  and  reliable  grid 
and  smoothing  criteria  are  unavailable  although  criteria  have 
been  developed  for  specific  algorithms  and  flow  conditions.  The 
present  effort  explores  criteria  that  may  be  more  generally 
useful  for  guiding  grid  adjustment  and  setting  smoothing  levels. 

Work  completed  during  an  earlier  contract  (titled  "Error  Norm 
Guided  Flow  Analysis" ,  contract  number  F49620-84C-0037 )  was 
summarized  in  a  symposium  paper  and  presentation  (Forester  and 
Strom,  1985).  These  results  dramatically  illustrated  the 
sensitivity  of  solution  results  to  grid  placement  and  density. 
Solutions  were  obtained  for  supersonic  flow  over  an  axisymmetric 
conical  afterbody  with  a  blunt  base  containing  a  centered 
propulsive  jet.  Comparisons  were  made  between  computed  and 
experimental  results  for  base  pressure,  separation  length, 
afterbody  pressure  distribution,  and  flow  field  structure.  The 
numerical  solutions  were  found  to  be  sensitive  to  the 
computational  grid  structure  and  the  turbulence  model.  Error 
norms  were  applied  to  aid  the  detection  of  inappropriate  grid 
choices.  The  best  results  were  obtained  with  adaptive  grids  that 
tracked  both  shear  layers  (i.e.,  due  to  the  internal  and  the 
external  flow)  and  a  turbulence  model  based  on  the  local  flow 
features.  The  results  of  this  work  suggested  that  artificial 
dissipation  as  defined  for  the  MacCormack  explicit  algorithm 
could  be  used  as  a  guide  to  grid  placement. 

A  variety  of  error  norms  were  explored  during  the  early  phases  of 
this  contract  (Forester,  1986,  2/87,  and  8/87).  This  early  work 
identified  ADR  (discussed  in  this  report)  as  a  viable  error 
monitor . 


2.  ERROR  SOURCES  AND  SMOOTHING 

When  a  computational  fluid  dynamic  (CFD)  analysis  is 
performed  the  objective  is  to  solve  a  set  of  partial  differential 
equations  ( PDE )  such  that 

PDE  -  0  (1) 
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Since  most  PDEs  are  too  complex  to  be  solved  exactly  it  is 
necessary  to  solve  an  approximation  to  the  PDE .  A  common 
approach  is  to  solve  a  set  of  finite  difference  equations  ( FDE ) 
such  that 

FDE  -  0  (2) 

but 

PDE  -  FDE  +  e  ( 3  ) 

Where  e  represents  nonphysical  source  terms  due  to  truncation, 
residual  and  roundoff  errors.  This  implies  that  the  equations 
actually  being  solved  are 

PDE  «  c  (4) 

Comparing  Equation  1  with  Equation  4  it  is  apparent  that  the 
resulting  solution  (obtained  by  solving  FDE  -  0)  contains  errors 
that  are  grid,  convergence  and  machine  dependent. 

The  flow  field  equations  and  their  boundary  conditions  are 
precise,  exact  and  continuous.  The  process  of  numerically 
simulating  these  equations  is  not  precise,  exact  or  continuous. 
Numerical  solutions  of  these  field  equations  involve  solving 
equations  which  have  additional  properties  that  are  not  contained 
by  or  related  to  the  continuous  equations.  These  additional 
properties  are  characteristic  of  the  algebraic  structure  (i.e., 
FDEs )  used  to  model  the  continuous  equations. 

These  additional  spurious  properties  have  been  catagorized  as 
dispersive  and  dissipative  behavior  in  the  solutions.  Examples 
of  spurious  behavior  are  nonphysical  mass  sources  and  sinks; 
negative  density  or  temperature;  and  the  decay  of  entropy  or 
increases  in  available  total  pressure.  Dispersive  errors  result 
in  nonphysical  oscillations  or  ringing  in  the  solution  to  the 
system  of  equations.  Dissipative  behavior  is  manifested  by  poor 
acuity  of  steep  gradients  or  by  damped  peak  amplitudes  of 
dependent  variables  or  their  derived  quantities.  The  combined 
effect  of  dispersive  and  dissipative  errors  is  often  called 
diffusion  error.  All  numerical  algorithms  whether  finite 
difference,  finite  volume  or  finite  element,  have  a  truncation 
error  and  thus  an  intrinsic  diffusion  unique  to  a  particular 
algorithm.  In  addition,  artificial  dissipation  is  often 
explicitly  added  to  an  algorithm  either  to  achieve  solution 
stability  or  to  minimize  undesirable  dispersive  errors  near 
regions  of  high  gradients  in  flow  properties.  The  combined 
effect  of  the  intrinsic  diffusion  and  the  explicitly  added 
artificial  dissipation  is  referred  to,  in  this  report,  as 
artificial  diffusion. 

If  the  numerical  solution  is  dispersive  near  a  high  gradient 
region  of  the  flow  the  effect  of  increasing  the  viscosity  is  to 
reduce  or  eliminate  (i.e.,  smooth)  the  local  oscillations  in  the 
solution.  Conversely,  reducing  the  viscosity  accentuates  the 
oscillations.  The  dispersive  numerical  error  is  counteracted  by 


the  viscous  dissipation  (i.e.,  has  an  opposite  effect  on  the  flow 
solution).  Since  the  source  term  on  the  right  hand  side  of 
Equation  4  typically  varies  (both  in  magnitude  and  sign)  from 
point-to-point  throughout  the  solution  domain,  judicious  addition 
of  numerical  viscosity  (or  dissipation)  in  the  vicinity  of  the 
high  gradient  region  will  reduce  or  eliminate  the  dispersive 
error.  In  effect  the  source  term  implied  by  Equation  4  is 
locally  reduced  or  eliminated. 

There  is  a  difference  in  the  magnitude  of  smoothing 
requirements  between  linear  and  nonlinear  regions  of  the  flow 
analysis  domain.  For  example,  regions  of  the  flow  featuring 
shocks  (nonlinear)  require  enormous  smoothing  levels.  Since 
shock-like  discontinuities  abound  in  large  numbers  during  the 
early  stages  of  the  convergence  process,  high  levels  of  smoothing 
are  required  in  most  of  the  flow  analysis  domain. 

As  the  converged  solution  is  approached,  the  smoothing  levels 
must  be  low  except  in  shock  regions  and  poorly  resolved  flow 
regions  that  behave  like  discontinuities.  These  regions  require 
locally  higher  smoothing  to  avoid  unacceptable  dispersive  errors. 
Locally  higher  smoothing  can  however  be  an  unacceptable 
dissipative  error.  Complexity  is  inherent  in  the  widely  varying 
demands  for  smoothing  during  the  convergence  process. 

Artificial  diffusion  can  smooth,  stabilize  and  enhance  the 
convergence  characteristics  of  NS  solution  algorithms  such  as  the 
MacCormack  predictor-corrector  scheme.  This  is  done  by 
constructing  a  smoothing  flux  that  is  added  explicitly  to  the  raw 
flux  in  the  predictor-corrector  cycles.  The  impact  on  the 
solution  accuracy  of  this  smoothing  operator,  relative  to  grid 
and  residual  effects,  has  not  been  well  defined  for  the  NS 
equations.  Model  equation  studies  from  the  follow-on 
investigation  suggest  that  improvements  in  solution  accuracy  can 
be  obtained  with  this  type  of  smoother. 

The  present  study  was  aimed  at  investigating  a  variety  of 
error  monitors  including  artificial  dissipation  to  direct  the 
analysis  process  for  achieving  "good"  grids  —  those  grids  giving 
solutions  suitable  for  engineering  application.  Care  was  taken 
to  insure  that  the  solutions  were  well  converged  and  that 
roundoff  error  was  negligible.  Therefore,  the  solution  errors 
discussed  in  this  report  are  dominated  by  the  truncation  (i.e, 
grid  dependent)  errors.  The  contract  and  follow-on  studies  were 
focused  on  truncation  error  and  smoothing  effects  (relative  to 
the  truncation  error). 


3.  ERROR  MONITORS 

During  the  contract  new  approaches  to  streamline  and  simplify 
the  process  of  selecting  an  analysis  grid  were  explored. 

Measures  of  solution  error  were  investigated,  including  total 
pressure  error  and  the  artificial  diffusion  ratio  defined  in  the 
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following  paragraphs.  Cross-correlation  of  these,  plus  other 
error  measures  such  as  local  and  global  conservation  error,  have 
been  examined  to  single  out  the  best  guides  for  selecting  grid 
and  smoothing  levels. 

Many  methods  have  been  used  to  select  the  local  level  of 
artificial  dissipation  or  smoothing  used  to  minimize  dispersive 
errors  in  a  flow  solution.  These  methods  range  in  complexity 
from  a  global  increase  in  the  bulk  viscosity  to  flux  corrected 
transport  (Oran  and  Boris,  1987).  The  MacCormack  (1975)  smoother 
(described  in  section  3.1  below)  senses  two  grid  interval 
solution  oscillations  (usually  a  dispersive  error  effect).  The 
artificial  dissipation  level  is  set  at  a  "high"  value  if  a 
dispersive  error  effect  is  present.  Because  the  MacCormack 
smoother  is  designed  to  sense  dispersive  error,  it  can  also  be 
used  as  a  guide  to  grid  adjustment. 

Artificial  smoothing  (MacCormack  method)  and  its  effect  on 
the  solution  accuracy  were  examined  in  this  study.  Truncation 
error  is  usually  the  dominant  error  source.  The  artificial 
diffusion  is  monitored  and  utilized  to  select  grid  densities  that 
ensure  accuracy  of  the  analysis  process.  It  is  shown  that  the 
artificial  diffusion  is  large  on  inadequate  or  coarse  grid 
solutions  where  flow  gradients  are  large.  This  information  is 
used  to  adjust  grid  density  until  the  artificial  diffusion  is 
minimized  and  an  accurate  solution  is  obtained. 

A  measure  of  the  artificial  diffusion  is  the  ratio  of  the 
explicit  smoothing  flux  (computed  with  a  MacCormack  smoother)  to 
the  total  flux  through  computational  cell  faces.  This  ratio  is 
called  the  artificial  diffusion  ratio  (ADR).  By  definition,  this 
ratio  should  be  insignificant  everywhere  in  the  flow  analysis 
domain  except  in  high  gradient  regions  where  the  mesh  is  too 
coarse.  For  example,  shock  waves  and  shear  layers  must  have  low 
levels  of  artificial  smoothing  or  small  values  of  ADR  relative  to 
the  peak  values  of  ADR.  Coarse  grid  simulations  lead  to  the 
treatment  of  each  of  these  continuous  regions  as  singularities. 

An  ADR  can  be  generated  for  each  dependent  variable  of  the 
compressible  time-averaged  Navier-Stokes  equations.  The 
dependent  variables  (in  2D)  are  two  components  of  momentum,  total 
energy  and  mass  density.  These  dependent  variables  result  in 
four  artificial  diffusion  ratios  —  ADRp ,  ADRq,  ADRh  and  ADRR, 
respectively.  A  composite  of  these  individual  error  monitors  may 
be  constructed  by  averaging  the  individual  error  monitors.  All 
error  measures  are  normalized  to  achieve  peak  values  not 
exceeding  unity. 

The  error  measures  that  have  been  explored  for  the  contract 
are  defined  in  the  next  section.  Finite  difference  expressions 
are  used  for  this  purpose. 
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3.1  Mathematical  Definition 


The  concept  of  error  analysis  for  numerical  approximations  of 
the  NS  equations  can  be  illustrated  using  a  model  equation  called 
the  Burgers  equation.  The  generalized  Burgers  equation  may  be 
written  in  a  form  consistent  with  the  usual  representation  of  the 
NS  equations  as  follows: 


3U  9E 
3t  3x 

where 


0 


(5) 


U-u,  E-F-yt/fp  F-cu  +  (6) 

u  is  the  velocity,  //  is  the  fluid  viscosity,  and  c  and  b  are  free 
parameters.  Note  that  if  b  -  0,  the  linearized  Burgers  equation 
results  and  if  c  -  0  and  b  -  1,  the  nonlinear  Burgers  equation  is 
obtained.  If  /j  -  0  the  equation  is  inviscid  and  if  b  -  0  and  u  ■ 
0  the  familiar  wave  equation  is  recovered. 

The  finite-difference  expression  for  MacCormack's  explicit 
predictor-corrector  method  with  a  "product"  smoother  is  given  by 
(Anderson  et  al.,  1984). 
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where  S  represents  the  artificial  dissipation  flux  at  a  point. 
For  the  NS  equations  S  is  defined  as 


sj  -  t  (|uj|  +  a?) 
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Note  that  when  solving  Burgers  equation,  the  pressure,  P,  is 
not  a  dependent  variable  and  the  speed  of  sound,  a,  is  not  a 
solution  parameter.  Therefore,  for  Burgers  equation,  velocity  is 
used  in  place  of  pressure  and  the  wave  speed,  c,  is  used  in  place 
of  "a"  in  the  definition  of  the  artificial  dissipation  flux. 
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An  additive  artificial  viscosity  is  thus  defined  as 


a 


.i  -  e  Ux)  r  K.l  -  2ui  +  ui-ll 


where 


r  - 


( |u. 1  +  c  ) 


(ui+l  +  2u.  +  ui_1 ) 


(11) 

(12) 


The  ax  is  included  in  the  definition  of  the  artificial  viscosity 
so  that  S  -  //a  3u/3x  and  is  comparable  to  the  viscous  term  v 
3u/3x  contained  in  E.  Also  note  that  the  sign  on  the  smoothing 
terms  in  Equations  7  and  8  have  been  selected  so  that  they  have 
the  same  sign  as  the  viscous  term  in  E.  This  smoother  is 
triggered  by  oscillations  in  a  numerical  solution  (a  third  order 
or  dispersive  effect).  A  monatonic  variation  of  u  over  the  i-1 
to  i+1  grid  interval  results  in  a  low  value  of  pa  A  two  grid 

interval  oscillation  in  u  results  in  a  high  value’of  ua  With 
an  appropriate  setting  for  e  the  dissipation  added  by  tfie 
smoother  will  damp  the  third  order  dispersion  error. 


The  artificial  diffusion  ratio,  ADR,  is  defined  as  follows: 


ADR  - 


+i 


<SU1 


(13) 


Note  that  the  magnitude  of  the  added  artificial  smoothing 
term  is  compared  with  the  3E/3x  term  such  that  ADR  must  range 
between  0  and  1.  This  definition  is  for  the  predictor  step.  A 
similar  definition  may  be  expressed  for  the  corrector  step,  but 
studies  have  shown  that  both  definitions  have  essentially  the 
same  characteristics.  ADR  based  on  the  predictor  step  was  used 
for  the  present  effort. 


The  ADR  may  be  defined  similarly  for  the  continuity  equation, 
the  two  components  of  the  momentum  equation,  and  for  the  energy 
equation.  These  quantities  are  called  ADRR,  ADRp,  ADRG,  and 
ADRe,  respectively. 


3.2  Minimizing  Errors 

The  equations  of  the  previous  section  constitute  a  system  of 
unknowns  whose  range  is  determined  by  the  number  and  distribution 
of  grid  points.  The  number  of  grid  points  is  chosen  rather 
arbitrarily  until  data  is  available  on  the  truncation  error. 

Once  this  data  exists,  a  systematic  procedure  for  minimizing  the 
error  using  grid  refinement  may  be  followed.  It  is  usually 
desirable  to  distribute  the  grid  so  that  the  truncation  error  is 
approximately  uniform  over  the  grid  (using  ADR  as  a  guide). 


8 


«• 


Once  the  grid  distribution  is  arrived  at,  the  accuracy  of  the 
numerical  approximation  can  be  increased  by  increasing  the  number 
of  grid  points.  If  this  is  done,  a  correlation  emerges  between 
the  truncation  error  and  ADR.  It  is  this  correlation  which  is 
critical  to  the  viability  of  ADR  as  a  useful  error  monitor. 


4.  TEST  PROBLEM 

To  test  the  error  monitors,  an  analysis  problem  was  required. 
It  was  desirable  to  select  this  problem  to  achieve  accurate 
computational  solutions  at  a  reasonable  cost.  High  quality  test 
data  was  also  needed  so  that  comparisons  could  be  made  between 
the  computed  and  measured  results.  The  inlet/aperture  flow  field 
test  problem  was  selected  because  the  test  data  was  of  high 
quality  and  because  the  geometry  was  two-dimensional.  Two- 
dimensional  NS  solutions  on  a  number  of  grids  were  feasible  for 
the  inlet  within  the  available  computer  budget. 


4.1  Geometry  and  Flow 

A  supersonic  inlet  throat  region  geometry  and  the  flow 
physics  associated  with  an  external  compression  inlet  with  throat 
bleed  flow  is  illustrated  in  Figure  1.  The  inlet/aperture 
approach  flow,  parallel  to  the  centerbody  ramp,  is  at  the  throat 
design  Mach  number  of  1.28.  The  bulk  of  the  flow  in  the 
streamtube  spanning  the  gap  between  the  cowl  lip  and  ramp  crown 


Some  spillage  flow  may  occur  at  the  cowl  lip,  depending  on 
the  engine  face  flow  and  the  bleed  slot  flow  requirements.  The 
cowl  generates  a  shock  wave  that  is  positioned  upstream  of  the 
lip  and  impinges  on  the  ramp  boundary  layer  downstream,  but  near 
the  crown  of  the  ramp.  The  ramp  shock-wave/boundary-layer 
interaction  is  stabilized  in  this  position  by  bleeding  flow  into 
a  slot  opening  in  the  ramp.  This  slot  is  downstream  of  the  crown 
of  the  ramp.  The  crown  of  the  ramp  generates  an  expansion  fan 
which  locally  accelerates  the  flow  to  a  higher  Mach  number  (= 

1.7)  before  the  flow  arrives  at  the  slot  opening.  These  flow 
physics  generate  a  complex  interaction  of  the  cowl  shock  wave, 
ramp  boundary  layer  and  the  ramp  crown  expansion  fan.  These 
interactions  produce  a  free  shear  layer  across  the  slot  opening 
which  has  strong  viscous  generated  transverse  velocity  gradients 
and  strong  local  inviscid  transverse  velocity  gradients  in  the 
supersonic  tongue.  Further,  this  flow  region  is  strongly 
influenced  by  the  slot  geometry  and  amount  of  bleed  flow.  These 
flow  and  geometry  features  generate  higher  gradients  in  flow 
properties  in  the  longitudinal  direction.  The  locations  of  these 
longitudinal  contortions  are  dictated  by  the  downstream  lip  of 
the  slot  opening.  This  lip  produces  a  strong  shock  wave  ( M^  = 
1.7)  in  the  slot  region.  It  is  a  stronger  shock  than  the  cowl 
shock  with  high  local  curvature.  This  curvature  adjusts  the  flow 
field  in  the  neighborhood  of  the  slot  opening  so  that  the 
resultant  flow  downstream  of  the  slot  opening  is  subsonic  to 
match  the  static  pressure  of  the  remaining  flow  captured  by  the 
inlet . 

The  slot  shock  and  the  shear  layer  interaction  are  sensitive 
to  the  supersonic  flow  originating  from  the  ramp  crown.  Flow 
solutions  near  the  slot  were  very  sensitive  to  solution  error. 
Incomplete  Prandtl-Meyer  expansion  in  the  numerical  solution 
dramatically  reduces  slot  shock  strength,  and  in  the  extreme 
case,  produces  no  slot  lip  shock.  Over-expansion  in  the  solution 
leads  to  misdirecting  the  slot  free-shear  layer  to  the  slot  lip. 
Either  of  these  mismatches  lead  to  incorrect  slot  pressure, 
either  too  high  or  too  low  respectively.  In  turn,  the  slot  flow 
entrance  separation  point  is  affected  dramatically  by  the 
pressure  in  the  bleed  flow  opening.  The  nature  of  the 
recirculation  flow  field  impinging  upon  the  bleed  entrance  free 
shear  layer  is  subject  to  change  according  to  the  separation 
point  which  induces  feedback  upon  the  bleed  entrance  pressure. 
Experimental  results  for  the  inlet  were  thought  to  be  sensitive 
to  experimentally  undefined  boundary  conditions. 

Turbulence  in  the  bleed  entrance  and  in  the  bleed  cavity 
affect  the  flow  pressure  at  the  bleed  entrance.  The  ratio  of  the 
length  of  the  bleed  flow  opening  and  the  thickness  of  the  free 
shear  layer,  plus  the  pressure  gradient  in  the  free  shear  layer, 
influence  the  structure  of  the  turbulence  in  the  bleed  cavity  and 
in  the  free  shear  layer.  Favorable  pressure  gradients  accelerate 
the  wall  boundary  layer  at  the  crown  of  the  ramp.  Mixing  zone 
intensity  is  excited  by  adverse  pressure  gradients  which  prevail 
between  the  bleed  cavity  and  the  supersonic  flow  regions. 


10 


Feedback  between  the  turbulence  and  the  pressure  field  can  cause 
unsteadiness  or  instability  of  the  flow  at  the  bleed  slot 
entrance.  The  intensity  of  the  mixing  can  be  damped  or  amplified 
by  the  unsteadiness.  The  mixing  process  may  include  at  least 
five  types  of  flow: 

1)  nearly  laminar  flow  of  the  wall  boundary  layer  at  the  ramp 
crown  and  downstream  of  the  separation  point, 

2)  transition  to  steady  flow  mixing  between  the  fully 
developed  supersonic  flow  and  the  bleed  cavity  flow  field, 

3)  possible  vigorous  flow  mixing  between  the  established 
bleed  flow  opening  free  shear  layer  and  the  bleed  cavity 
flow  field, 

4)  possible  resonance  phenomena  of  the  free  shear  layer  and 
the  bleed  cavity  flow  field,  and 

5)  fully  developed  mixing  in  boundary  and  free  shear  layers. 

The  occurrence  of  these  types  of  flow  phenomena  depends  upon 
flow  transition  phenomena.  Model  scale  inlets  and  full  scale 
inlets  can  have  large  differences  in  the  location  and  streamwise 
extent  of  transition.  For  the  model  scale  inlet,  examined  in  the 
present  study,  only  flow  types  1,  2  and  5  have  been  identified. 
Numerical  experiments  indicate  that  only  type  5  flow  dominates 
the  wall  boundary  layers  and  free  shear  layers. 


4.2  Computational  Approach 

The  computational  approach  in  the  present  study  is  an 
extension  of  the  work  of  Peery  and  Forester.  It  uses  a 
conservation-based  body-fitted  adaptive  grid  model  of  the  thin- 
shear-layer  formulation  of  the  compressible,  Reynolds-averaged, 
Navier-Stokes  equations  together  with  mass  and  energy 
conservation  equations. 

Control  of  the  residual  errors  is  achieved,  by  an  artificial 
time  relaxation  approach  with  a  constant  CFL  criteria.  Steady 
state  is  achieved  by  asymptotic  time  relaxation.  Truncation 
errors  are  reduced  through  the  use  of  solutions  on  varying  grid 
densities  and  varying  grid  distributions.  The  formulation  of  (a) 
the  governing  thin-shear-layer  equations,  (b)  the  finite  volume 
explicit  predictor-corrector  finite  difference  algorithm,  (c)  the 
boundary  conditions,  (d)  the  two-layer  algebraic  turbulence  model 
with  its  associated  wall  functions,  and  (e)  the  mesh  generator 
and  the  adaptive  mesh  mover  are  detailed  by  Campbell  and  Forester 
(1985)  and  the  associated  references.  The  procedure  for  design 
application  of  this  code  is  given  by  Campbell  et  al.  (1984). 

The  computer  program  allows  three  coupled  computational 
regions.  In  the  present  study,  computational  region  one  (see 
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Figure  2)  is  the  slot  cavity  flow  region,  region  two  is  the 

ramp/cowl/engine-face  flow  region,  and  region  three  is  the  region 

above  the  cowl.  The  grid  blocks  have  point  and  slope  continuity  I 

at  each  interface.  This  simplifies  the  boundary  condition 

treatment  needed  at  the  surface  of  grid  blocks.  Freestream 

boundary  conditions  are  specified  on  the  left  side  of  regions  two 

and  three.  The  grid  is  body-fitted  to  the  ramp,  slot,  and  cowl. 

As  shown  in  Figure  2,  the  length  scales  of  the  grid  intervals 

vary  widely  over  the  analysis  domain.  The  smallest  grid  j 

intervals  are  generally  located  in  critical  regions  as  follows: 

o  boundary  layers, 

o  free  shear  layers, 

i 

o  rapid  compression/expansion  regions, 
o  shocks,  and 
o  stagnation  regions. 

To  simplify  the  notation  for  grid  size  definition,  grid  sizes 
are  defined  by  a  cluster  of  numbers  separated  by  commas.  The 
numbers  between  the  commas  are  the  interval  counts  in  x  by  y 
directions  of  the  grid  in  region  1,  region  2  and  region  3, 
respectively.  The  grid  sizes  employed  are  (10x5,  27x17,  22x12), 

(21x10,  54x34,  44x24),  and  (42x18,  106x66,  46x26).  These  three 
grid  sizes  are  labeled  coarse,  medium  and  fine,  respectively. 

The  influence  of  grid  on  resolution  is  considered  relative  to 
accuracy  produced  by  pairs  of  grids  (coarse/medium,  coarse/fine, 
medium/fine.)  Only  the  coarse  and  fine  grid  results  are  shown. 


Analysis  Domain 


Figure  2.  Fine  Grid 
-  12  - 


Slot  Region 


4.3  Computational  Results 


Many  parameters  could  be  studied  to  explore  error  norm 
guides.  The  parameters  selected  for  this  study  were  the  total 
pressure,  artificial  diffusion  ratio,  and  Mach  number.  All  of 
these  parameters  except  for  the  artificial  diffusion  ratio  are 
traditional  parameters  for  error  assessment.  The  artificial 
diffusion  ratio  is  shown  in  conjunction  with  the  Mach  contours 
and  does  not  replace  these  because  physical  features  simply  are 
not  revealed  by  the  artificial  diffusion  ratio.  In  fact,  the 
artificial  diffusion  ratio  should  be  void  of  physical  features 
except  near  shocks.  Mach  contours  show  physical  features 
including  shocks. 

Figures  3  through  7  relate  solution  accuracy  to  smoothing 
coefficient  level  and  to  grid  placement  through  test  data 
comparisons.  Figure  1  shows  the  traverse  station  for  the  total 
pressure  ratio,  Pt/Pt«,  in  the  throat  region  of  the  inlet. 

Figure  3  shows  the  improvement  in  the  total  pressure  with  respec 
to  grid  refinement  and  with  respect  to  reducing  the  artificial 
smoothing. 


See  Figure  1 .  for  Transverse  Station 


Figure  3.  Total  Pressure  Ratio  in  Throat  Region 


if 


Note  that  the  degradation  of  the  total  pressure  results  from 
excessive  smoothing  or  from  too  coarse  a  grid.  Also  note  that 
the  results  improve  with  respect  to  grid  refinement  even  with 
abnormally  high  smoothing  coefficients.  Ultimately  grid 
refinement  is  the  critical  issue.  However,  the  efficiency  of 
using  a  particular  grid  is  improved  by  setting  smoothing 
coefficients  near  the  stability  limit,  rather  than  to  maximize 
convergence.  Further  examples  of  this  behavior  are  now 
discussed. 

Figures  4  through  7  are  comparisons  of  solution  results  for 
selected  grids  and  smoothing  coefficients.  Figures  4  and  5  show 
the  effect  of  smoothing  level  on  ADR  for  coarse,  and  fine  grids, 
respectively.  Figures  6,  7,  and  8  are  Mach  number  contours. 
Figures  4  and  5,  and  Figures  6  and  7  show  the  effect  of  grid 
density.  Comparison  of  Figures  4  and  5  with  Figures  6  and  7  show 
that  ADR  rises  sharply  with  increased  smoothing  levels  and  with 
increased  grid  coarseness.  Note  that  when  ADR  is  above  0.01 
(except  for  shocks  where  ADR  should  be  about  0.01),  too  high  a 
smoothing  level  or  too  coarse  a  grid  is  indicated.  In  these 
regions,  the  grid  must  be  refined  or  the  smoothing  level  must  be 
reduced. 

It  is  possible  to  generate  a  composite  effect  of  all  of  these 
error  sources  on  a  particular  grid.  The  grid  used  for  this 
purpose  is  shown  in  Figure  2.  Figure  8  shows  an  example  of  Mach 
number  contours  for  the  aperture  region.  Note  the  agreement 
between  the  shadowgraph  of  the  flow  field  for  an  experimental 
test  of  this  inlet  and  the  predicted  result  for  the  same  flow 
field  ( Figure  8 ) . 


Figure  4.  Coarse  Grid  Artifical  Diffusion  Ratio  Contours 
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Shadow  Graph 


Predicted  Mach  Contours 


Figure  8.  Comparison  of  Experiment  with  Analysis 
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5.  RELATED  RESEARCH 


#• 


Because  of  the  promising  nature  of  the  contract  results,  a 
follow-on  study  was  conducted  to  address  several  remaining 
questions.  The  follow-on  study  was  performed  with  BAS  IRAD 
funds.  In  particular,  the  questions  addressed  were  as  follows: 

o  Does  the  directionality  of  ADR  (i.e.,  the  u  and  v  momentum 
ADR,  ADRp  and  ADRq)  provide  guidance  to  the  direction  of 
the  required  grid  refinement? 

o  Since  the  previous  results  indicate  some  error  still  exists 
between  the  predicted  and  measured  total  pressures,  can 
grid  refinements,  that  achieve  the  suggested  ADR  levels, 
reduce  this  error? 

o  Is  ADR  useful  for  algorithms  other  than  the  MacCormack 

explicit  algorithms  studied  thus  far?  Specifically,  is  ADR 
useful  for  the  Beam-Warming  algorithm  used  in  the  popular 
ARC  NS  codes  (Pulliam,  1984)  and  the  PARC  NS  codes  (Cooper, 
1987)? 

o  What  is  the  relationship  between  ADR  magnitude  and  the 
solution  error  magnitude? 

o  Can  ADR  be  used  with  algorithms  that  use  intrinsically 
added  smoothing  such  as  the  MacCormack  implicit  algorithm 
(MacCormack,  1985)? 

The  usefulness  of  ADR  to  guide  grid  adjustment  in  a  grid 
direction  (directionality)  and  the  question  of  whether  grid 
adjustment  for  the  inlet  problem  to  reduce  local  ADR  levels  would 
improve  agreement  between  the  computed  and  measured  flow 
properties  were  addressed  by  continuing  the  contract  study 
approach  -  application  of  the  NS  analysis  to  the  inlet  test 
problem.  Work  and  results  from  this  are  reported  in  Section  5.1 
below.  The  usefulness  of  ADR  for  other  algorithms  and  the 
relationship  between  ADR  and  solution  accuracy,  were  addressed 
through  model  equation  studies  using  Burgers  equation  as  an 
analog  to  the  NS  equations.  Work  and  results  from  this  are 
reported  in  Section  5.2. 


5.1  ADR  Directionality  and  Grid  Refinement 

Additional  grid  refinement  work  was  performed  for  the  test 
case,  discussed  in  Section  4.3,  using  BAS  IRAD  funds  (Baltar, 
1988).  Details  of  this  study  are  included  in  Appendix  A  and  are 
summarized  in  this  section.  The  ADR  measures  were  used  to  guide 
this  grid  refinement  study. 
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During  the  contract  effort  ADR  was  identified  as  a  useful 
parameter  for  guiding  grid  refinement  and  establishing  solution 
accuracy.  The  objective  of  the  IRAD  effort  summarized  in 
Appendix  A  was  to  investigate  if  ADR  is  a  useful  guide  for 
localized  grid  refinement  of  a  complex  flow  field  analysis. 
Improvement  in  flow  prediction  was  evaluated  by  comparing  the 
calculated  total  pressure  with  experimentally  measured  total 
pressure  in  the  throat  region.  Regions  of  high  ADR  were  shown  to 
be  reduced  in  size  when  more  mesh  was  added  in  the  direction 
normal  to  the  gradient  or  in  a  way  to  reduce  the  cell  aspect 
ratio  of  the  mesh.  The  agreement  between  computed  and 
experimentally  measured  total  pressure  profiles  improved  about  2% 
when  the  recommended  levels  of  ADR  were  achieved  throughout  the 
computed  flow  field.  The  remaining  disagreement  could  be  caused 
by  improper  ADR  thresholds,  improper  geometry  resolution,  or  by 
experimental  uncertainties.  Preliminary  results  from 
investigation  into  the  use  of  ADR  to  determine  the  direction  of 
the  required  grid  refinement  were  promising  but  additional  study 
is  needed. 

During  the  course  of  this  study  a  valuable  graphical 
technique  for  using  ADR  to  guide  grid  refinement  was  developed. 
Figure  9  is  a  representative  ADRE  contour  plot.  Note  the  three 
color  bands.  Since  ADR  thresholds  of  0.01  and  0.001  were  being 
evaluated  the  color  bands  represent  regions  where  the  ADR  is 
greater  than  0.01  (the  red  band),  between  0.01  and  0.001  (the 
yellow  band),  and  less  than  0.001  (the  green  band).  The 
procedure  for  using  these  plots  is  to  refine  the  grid  in  the  red 
and  yellow  regions  first.  The  green  regions  were  assumed  to  have 
an  adequate  grid  density.  Once  the  majority  of  the  flow  domain 
had  relatively  low  ADR  values  (illustrated  by  the  large  green 
region  in  Figure  9)  the  agreement  with  the  data  was  good.  As 
shown  in  Figure  10,  the  agreement  between  the  predicted  and 
measured  total  pressures  in  the  throat  region  are  within  2%  and 
for  much  of  the  traverse  the  agreement  is  nearly  perfect.  Also 
note  the  crispness  of  the  shock  definitions  in  the  Mach  contour 
plot . 

5.2  Model  Equation  Studies 

Two  problems  were  evident  with  continuation  of  the  contract 
approach  addressing  the  usefulness  of  ADR  for  other  algorithms 
and  establishing  the  relationship  between  solution  error  and  ADR. 
First,  exploring  the  usefulness  of  ADR  with  other  algorithms  was 
too  complex  to  be  practical  (in  the  context  of  the  short  term 
follow-on  study)  Uoing  solutions  of  the  NS  equations.  Second, 
establishing  the  level  of  accuracy  of  solutions  to  the  NS 
equations  is  difficult  because  analytic  solutions  do  not  exist 
for  cases  of  interest.  Use  of  experimental  results  to  represent 
an  analytic  (exact)  solution  to  the  NS  equations  is  dangerous 
because  of  experimental  uncertainties. 
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Figure  10.  Inlet-Aperture  Analysis:  Mach  Number  Contours  and  Total  Pressure  Profile 


The  approach  selected  for  the  follow-on  study  was  as  follows: 

o  Use  Burgers  equation  as  a  simple  analog  of  the  NS 
equations . 

o  Solve  Burgers  equation  using  the  MacCormack  and  the  Beam- 
Warming  algorithms  both  with  and  without  artificial 
dissipation  modeling. 

o  Establish  the  accuracy  levels  of  the  numerical  solutions 
through  comparisons  between  analytic  and  numerical  results. 

o  Complete  modified  equation  analyses  of  the  MacCormack  and 
Beam-Warming  FDE  approximations  to  Burgers  equation  to 
address  why  ADR  might  be  a  measure  of  solution  accuracy. 

The  advantages  of  this  approach  are  that  Burgers  equation  is 
recognized  as  an  analog  to  the  NS  equations,  analytic  solutions 
are  available  to  Burgers  equation  so  that  the  accuracy  of  the 
numerical  solutions  can  be  established  without  ambiguity,  the 
MacCormack  and  Beam-Warming  algorithms  are  currently  the  most 
popular  explicit  and  implicit  methods,  and  programming  the  two 
methods  is  simple  enough  to  be  practical  within  the  timeframe  of 
a  short  study.  Further,  while  a  successful  result  with  ADR  in 
the  follow-on  study  doesn't  guarantee  the  success  of  using  ADR 
with  the  NS  equations,  failure  would  make  success  with  the  NS 
equations  improbable.  Success  here  is  defined  as  establishing  a 
relationship  between  ADR  and  solution  error  and  demonstrating  the 
usefulness  of  ADR  for  the  Beam-Warming  and  MacCormack  algorithms. 


5.2.1  Model  Equation  Analysis 

Using  Burgers  equation  as  an  analog  to  the  NS  equations, 
modified  equations  were  developed  for  the  linear  viscous  Burgers 
equation  for  both  the  MacCormack  explicit  method  (used  for  the 
present  contract  work)  and  the  Beam-Warming  implicit  method  (used 
in  other  popular  NS  codes  such  as  ARC  and  PARC).  The  details  of 
this  portion  of  the  follow-on  IRAD  study  (Paynter,  1988)  are 
included  in  Appendix  B  and  are  summarized  below. 

The  modified  equation  is  the  finite  difference  expression 
written  as  a  sum  of  the  partial  differential  equation  and  the 
truncation  error  with  the  lead  time  derivative  terms  in  the 
truncation  error  replaced  by  spatial  derivatives  through  a 
manipulation  of  this  equation.  Modified  equations  were  developed 
for  the  linear  viscous  Burgers  equation  for  both  the  MacCormack 
explicit  and  the  Beam-Warming  implicit  algorithms.  It  was  found 
that  the  lead  truncation  error  term  with  the  MacCormack  algorithm 
is  dispersive  and  that  the  lead  term  with  the  Beam-Warming 
algorithm  is  dissipative. 
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The  artificial  diffusion  is  in  effect  a  source  term  added  to 
the  PDE  being  solved.  The  magnitude  of  the  source  term  is  at 
least  a  qualitative  measure  of  the  dispersive  solution  error. 
Because  of  the  non-linear  nature  of  the  flow  equations,  a 
quantitative  relationship  between  solution  error  and  ADR  was  not 
established.  A  quantitative  relationship  between  solution  error 
and  ADR  would  be  algorithm  and  problem  dependent  due  to  the 
differences  in  the  lead  terms  of  the  truncation  error  from 
algorithm  to  algorithm.  This  is  supported  by  the  computed 
results  of  Section  5.2.2.  The  artificial  diffusion  should  be 
less  useful  for  algorithms  with  high  intrinsic  dissipation  since 
ADR  is  insensitive  to  intrinsic  dissipation  effects.  It  seems 
clear  that  a  quantitative  relationship  between  artificial 
diffusion  and  solution  accuracy  would  be  algorithm  and  problem 
dependent . 


5.2.2  ADR  as  a  Measure  of  Error  Magnitude 

Comparisons  between  analytic  and  CFD  solutions  of  Burgers 
equation  for  a  range  of  qrids  for  the  Beam-Warming  and  MacCormack 
algorithms  were  studied  as  part  of  the  follow-on  IRAD  study 
(Mayer,  1988).  The  details  of  this  are  reported  in  Appendix  C 
and  the  results  summarized  below. 

o  solution  accuracy  is  not  simply  related  to  ADR  thresholds, 

o  ADR  is  useful  for  defining  grid  regions  that  require 
refinement,  and 

o  ADR  is  less  useful  in  defining  regions  of  high  grid  error 
for  the  Briley-McDonald  algorithm  (which  is  equivalent  to 
the  Beam-Warming  algorithm  when  applied  to  Burgers 
equation)  than  for  the  MacCormack  explicit  algorithm. 

The  computer  program  developed  for  this  follow-on  study 
(Appendix  D)  could  be  readily  adapted  for  other  algorithms, 
including  those  with  intrinsic  dissipation.  Modified  equation 
analyses  could  be  performed  to  define  the  lead  truncation  error 
terras  for  other  algorithms  and  model  equations.  Analytic 
solutions  could  be  used  to  directly  evaluate  these  error  terms. 

A  detailed  study  of  these  error  terms  relative  to  the  definition 
of  ADR  should  provide  valuable  insight  and  guidance  to  using  ADR 
for  assessing  solution  accuracy. 


6.  RESEARCH  STATUS 


Several  significant  accomplishments  and  conclusions  have 
resulted  from  the  contract  effort  and  associated  Boeing  Advanced 
Systems  (BAS)  independent  research  and  development  (IRAD)  funded 
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research.  In  addition,  the  promising  results  from  the  current 
work  suggests  additional  research  to  develop  the  error  monitor, 
the  artificial  diffusion  ratio  (ADR),  into  a  tool  useful  for 
general  purpose  CFD  applications.  These  accomplishments, 
conclusions  and  recommendations  are  discussed  in  the  following 
sections . 


6.1  Significant  Accomplishments 

The  solutions  of  the  flow  field  in  the  aperture  region  of  an 
external  compression  inlet  with  bleed  and  spillage  flow  by 
Navier-Stokes  analysis  were  obtained  on  a  variety  of  grids  and  a 
new  approach  was  explored  to  guide  grid  adjustment.  Measures  of 
numerical  errors  in  the  analysis  process  were  explored  including 
ADR  for  mass,  energy,  and  momentum.  Correlation  of  these  error 
measures  show  that  ADR  provides  guidance  for  grid  and  smoothing 
level  selection.  The  application  of  ADR  leads  to  a  grid  choice 
that  yields  an  adequate  solution  to  the  flow  field.  Comparison 
of  this  solution  with  experimental  data  shows  good  agreement. 

The  promising  contract  results  led  to  a  follow-on  IRAD  study  to 
further  explore  the  usefulness  of  ADR  as  a  guide  to  grid 
adjustment . 

In  the  follow-on  study,  the  usefulness  of  ADR  was  further 
explored  for  the  inlet  test  problem  on  a  variety  of  grids.  The 
applicability  of  ADR  to  other  algorithms  and  the  relationship 
between  ADR  and  solution  accuracy  were  explored  in  model  equation 
studies.  The  inlet  studies  resulted  in  a  new  procedure  for 
graphical  display  of  the  ADR  data  and  some  work  toward 
establishing  the  usefulness  of  ADR  to  guide  grid  adjustment  in  a 
given  direction. 

Modified  equations  were  developed  for  the  linear  viscous 
Burgers  equation  for  both  the  MacCormack  explicit  and  the  Beam- 
Warming  implicit  algorithms.  This  development  clearly 
illustrates  the  nature  of  the  lead  truncation  error  terms  for 
these  algorithms.  The  nature  of  these  terms  indicate  how  ADR 
detects  solution  errors  for  MacCormack's  explicit  method.  In 
addition,  the  mathematical  development  of  these  terms  makes  it 
feasible  to  directly  evaluate  the  magnitude  of  these  errors  for 
test  problems  that  have  analytic  solutions. 

Analytic  and  CFD  solutions  on  a  range  of  grids  were  obtained 
for  the  viscous  Burgers  equations  for  both  the  Beam-Warming  and 
MacCormack  implicit  algorithms.  Comparisons  between  the  analytic 
and  CFD  results  were  useful  in  establishing  the  relationship 
between  solution  accuracy  and  ADR  for  a  given  algorithm.  The 
computer  program  written  for  the  model  equation  studies  could  be 
easily  extended  for  the  study  of  other  algorithms. 
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6.2  Conclusions 


The  major  conclusions  drawn  from  the  research  are  as  follows: 

o  The  use  of  ADR  to  pinpoint  grid  problems  has  been 

demonstrated.  Regions  of  the  flow  where  the  grid  is  too 
coarse  to  accurately  capture  gradients  results  in  noise  or 
ringing  in  the  solution.  This  noise  is  dispersive  in 
nature,  and  when  using  a  MacCormack  smoother,  results  in 
correspondingly  high  values  of  ADR.  Therefore,  regions  of 
the  flow  that  have  relatively  high  levels  of  ADR  are 
candidates  for  further  grid  refinement. 

o  Solution  accuracy  is  improved  and  the  high  ADR  values  are 
reduced  by  performing  localized  grid  refinement  based  on 
the  location  of  high  ADR  values. 

o  ADR  thresholds  as  direct  measures  of  solution  accuracy  have 
not  been  established.  There  does  exist  a  correlation 
between  ADR  and  solution  accuracy,  but  the  correlation  is 
more  complex  than  a  simple  one-to-one  relationship. 
Achieving  a  given  level  of  ADR  does  not  guarantee  a  given 
accuracy  level.  This  relationship  is  thought  to  be 
algorithm  and  problem  dependent. 

o  ADR  does  not  pinpoint  regions  of  grid  error  as  well  for  the 
Briley-McDonald  method  (comparable  to  the  Beam-Warming 
method  for  the  ID  Burgers  equation).  The  regions  of  high 
ADR  cover  a  greater  extent  than  the  regions  of  relatively 
high  solution  error.  This  means  that  ADR  for  the  Briley- 
McDonald  method  exaggerates  the  extent  of  regions  that 
require  grid  refinement. 

o  The  lead  truncation  error  term  is  dispersive  for  the 
MacCormack  algorithm  and  is  dissipative  for  the  Beam¬ 
warming  algorithm.  In  addition,  an  intrinsic  third  order 
dispersive  error  term  exists  for  the  Beam-Warming 
algorithm. 

o  The  artificial  dissipation  acts  as  a  source  term  added  to 
the  partial  differential  equation  being  solved.  For 
algorithms  with  high  dispersive  error  the  artificial 
dissipation  is  a  qualitative  measure  of  solution  accuracy. 

o  Artificial  dissipation  is  less  useful  for  algorithms  with 
high  intrinsic  dissipation. 

o  A  quantitative  relationship  between  artificial  dissipation 
and  solution  accuracy  is  algorithm  and  problem  dependent. 
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6.3  Recommendations 

Further  work  is  recommended  for  demonstrating  the  utility  of 
ADR  as  follows: 


o  The  results  of  the  modified  equation  analysis  of  the 
MacCormack  explicit  and  the  Beam-Warming  implicit 
algorithms  (Paynter,  1988,  Appendix  B)  makes  it  feasible  to 
directly  evaluate  the  lead  truncation  error  terms  for  an 
analytic  solution.  A  detailed  study  of  these  terms 
relative  to  the  present  (or  alternate)  definitions  of  ADR 
should  provide  valuable  insight  and  guidance  to  using  ADR 
for  assessing  the  solution  accuracy  of  these  algorithms. 

o  The  extension  to  solution  algorithms  that  use  intrinsic 
smoothing  (i.e.,  not  added  explicitly)  has  not  been 
demonstrated.  The  computer  program  written  for  the  present 
study  (Appendix  D)  may  be  readily  adapted  to  explore  the 
use  of  ADR  for  these  and  other  algorithms. 

o  Additional  test  problems  should  be  investigated  to  evaluate 
the  use  of  ADR  (for  algorithms  studied  with  the  above 
methods)  for  flows  with  different  flow  field 
characteristics  (i.e.,  flow  separation,  without  shocks, 
inviscid,  laminar,  etc.).  These  test  problems  should  be 
performed  for  configurations  with  high  quality  experimental 
data  suitable  for  CFD  validation. 

o  The  utility  of  ADR  should  be  tested  for  other  design 

analyses  of  aircraft  components  (both  model  scale  and  full 
scale ) . 
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APPENDIX  A: 

Artificial  Diffusion  Ratio  (ADR)  Guided  Grid  Refinement  Study 

by 

JY  Baltar 


1 .  Summary 

Forester  and  Tjonneland  (1988)  identified  the  artificial  diffusion  ratio 
(ADR)  as  a  useful  parameter  for  guiding  grid  refinement  and  establishing 
solution  accuracy.  Global  grid  refinement  and  artificial  viscosity  levels 
were  studied  in  their  paper.  The  objective  of  this  appendix  is  to 
investigate  if  ADR  is  a  useful  guide  for  localized  grid  refinement  of  a 
complex  flow  field  analysis.  Improvement  in  flow  prediction  is  evaluated  by 
comparing  the  calculated  total  pressure  with  experimentally  measured  total 
pressure  in  the  throat  region.  Regions  of  high  ADR  are  shown  to  be  reduced 
in  size  when  more  mesh  is  added  in  the  direction  normal  to  the  gradient  or  in 
a  way  to  reduce  the  cell  aspect  ratio  of  the  mesh.  The  agreement  between 
computed  and  experimentally  measured  total  pressure  profiles  improved  about 
2%  when  the  recommended  levels  of  ADR  were  achieved  throughout  the  computed 
flow  field.  The  disagreement  that  remains  may  be  caused  by  improper  ADR 
thresholds,  improper  geometry  resolution,  or  by  experimental  uncertainties. 
Preliminary  investigation  into  the  use  of  ADR  to  determine  the  direction  of 
the  required  grid  refinement  looks  promising  but  more  study  is  needed. 

«• 

2.  Introduction 

This  appendix  documents  the  work  performed  using  ADR  to  guide  a  grid 
refinement  study  on  the  inlet-aperture  problem  analyzed  by  Campbell  and 
Forester  (1985).  The  previous  work  (performed  by  PS  Hertel  of  Boeing 
Computer  Services  and  CK  Forester)  studied  the  changes  in  ADR  during  global 
grid  refinement  and  artificial  viscosity  level  studies  on  this  inlet  flow 
field.  The  results  were  reported  by  Forester  and  Tjonneland  (1988). 

The  geometry  and  flow  field  features  for  the  inlet-aperture  case  analyzed 
are  shown  schematically  in  Figure  A-l.  The  approach  flow  parallel  to  the 
ramp  surface  is  at  a  Mach  number  of  1.28.  A  cowl  shock  wave  is  established 
which  impinges  on  the  ramp  boundary  layer  near  the  crown  of  the  ramp.  The 
ramp  boundary-layer/shock-wave  interaction  is  stabilized  in  this  position  by 
bleeding  flow  into  a  slot  opening  in  the  ramp.  The  crown  of  the  ramp 
generates  an  expansion  fan  which  locally  accelerates  the  flow  before  the  flow 
arrives  at  the  slot  opening.  The  downstream  lip  of  the  bleed  slot  opening 
produces  a  strong  shock  wave  so  that  the  resultant  flow  downstream  of  the 
slot  is  subsor  c  to  match  the  static  pressure  of  the  remaining  flow  captured 
by  the  inlet.  This  highly  complex  flow  region  with  normal  shocks,  rapid 
expansions,  and  strong  shear  layers  presents  a  challenging  problem  for  CFD 
analysis. 
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A  new  approach  was  introduced  by  Forester  and  Tjonneland  (1988)  to  aid  in 
guiding  grid  refinement  studies.  The  parameter  used  to  guide  the  grid 
refinement  is  ADR.  ADR  is  defined  as  the  ratio  of  the  artificial  smoothing 
flux  to  the  total  flux  through  computational  cell  faces.  Values  of  ADR  of 
0.01  or  less  in  shock  waves  and  values  of  ADR  of  0.001  or  less  for  smooth 
flow  regions  were  recommended  by  Forester  (1987).  The  goals  of  this  study 
were  to  1)  achieve  the  recommended  levels  by  locally  adding  grid  points  where 
the  ADR  was  high  to  see  if  the  computed  flow  solution  could  be  brought  into 
agreement  with  experiment;  and  2)  to  determine  if  the  directionality  of  ADR 
could  be  used  to  guide  the  x-  versus  y-direction  grid  refinement. 
Improvement  in  flow  prediction  was  evaluated  by  comparing  the  calculated 
total  pressure  with  experimentally  measured  total  pressure  in  the  throat 
region.  Since  experimental  uncertainties  exist,  perfect  agreement  with  the 
data  was  not  expected.  The  data  does  provide  valuable  indications  of  grid 
related  errors. 

As  described  by  Forester  and  Tjonneland  (1988),  ADR  can  be  generated  for 
each  dependent  variable  of  the  flow  equations  that  are  solved  by  the  CFD 
code.  For  this  study,  the  ADR  for  the  energy  equation  and  the  two  momentum 
equations  were  calculated  and  analyzed  to  determine  if  they  could  be  used  to 
determine  where  and  in  which  direction  the  grid  needed  refinement.  These 
ratios  were  labeled  ADRg  (energy),  ADRp  (streamwise  momentum)  and  ADRq 
(normal  momentum)  after  the  standard  names  for  the  flux  vectors  in  the  2D 
Navier-Stokes  equations. 


3.  Results  and  Discussion 

Figure  A-2  shows  a  color  contour  plot  of  ADR  for  the  energy  equation 
(ADRg)  with  the  grid  overlayed.  This  data  is  from  the  final  fine  grid/low 
smoothing  result  shown  by  Forester  and  Tjonneland  (1988).  This  case  will  be 
called  the  Reference  Case  for  the  remainder  of  this  discussion.  Table  A-l 
gives  some  pertinent  information  about  this  case  and  about  the  other  cases 
discussed  in  this  appendix.  In  Figure  A-2,  the  areas  in  which  ADRp  is  above 
the  threshold  of  0.01  are  colored  red  while  the  regions  in  which  ADRg  is 
between  0.001  and  0.01  are  colored  yellow.  The  regions  of  the  flow  field 
which  are  below  the  threshold  of  0.001  are  colored  green.  This  color  scheme 
is  used  throughout  this  appendix  since  both  ADR  thresholds  may  be  shown  on 
one  plot.  The  ADRp  and  ADRq  color  contour  plots  for  this  case  are  shown  in 
Figures  A-3  and  A-4. 

Figure  A-5  shows  a  plot  of  the  total  pressure  profile  from  this  CFD 
solution  compared  with  the  experimentally  measured  values.  The  total 
pressure  measuring  plane  is  about  0.012  meters  or  12  nodes  downstream  of  the 
bleed  slot  lip  and  is  perpendicular  to  the  wall  of  the  constant  area  section 
of  the  duct  (see  Figure  A-l).  The  agreement  between  the  CFD  prediction  and 
experiment  is  fair  considering  the  complexity  of  the  flow  field.  However, 
several  regions  of  the  flow  indicate  a  need  for  improvement  to  accurately 
predict  the  local  total  pressure.  Notably,  the  region  between  0.0075  and 
0.03  meters  from  the  lower  wall  (which  is  strongly  influenced  by  the 
expansion-fan/shock-wave  interaction)  shows  a  disagreement  of  as  much  as  4%. 


During  this  study,  several  refinements  were  made  to  the  grid  that  is 
shown  in  Figure  A-2.  Since  the  0.01  threshold  regions  for  the  Reference  Case 
(Figures  A-3  and  A-4)  were  larger  for  ADRq  than  for  ADRp,  the  initial  grid 
refinements  were  made  by  adding  more  mesh  in  the  y-direction  (normal  to  the 
flow).  This  was  accomplished  while  holding  the  near  wall  cell  height 
constant  so  the  maximum  cell  height  at  the  center  of  the  duct  was  reduced 
significantly.  For  this  grid  refinement,  24  cells  were  added  in  the  y- 
direction  in  two  steps  for  a  total  of  90.  The  results  for  this  case  (Case 
la)  are  shown  in  Figures  A-6,  A-7  and  A-8.  Figures  A-6  and  A-7  show  that  the 
0.01  threshold  regions  in  the  cowl  and  slot  shocks  were  almost  removed  and 
the  0.001  threshold  regions  in  the  smooth  flow  regions  (particularly  in  the 
expansion  from  the  upstream  end  of  the  bleed  slot)  were  reduced 
significantly.  Figure  A-8  shows  that  the  flow  field  total  pressure  profile 
comparison  is  improved  so  that  the  maximum  disagreement  is  now  less  than 
about  2%.  This  improvement  can  be  attributed  to  the  increased  resolution  in 
the  region  of  the  expansion  and  the  shock  wave/expansion  interaction.  The 
increased  resolution  in  this  region  is  primarily  due  to  the  fact  that  the 
cell  aspect  ratio  was  reduced  to  near  unity  in  a  region  with  strong  gradients 
diagonal  to  the  mesh. 

A  further  refinement  in  the  y-grid,  called  Case  lb,  brought  the  total 
number  of  cells  in  the  y-direction  to  132  which  is  double  the  amount  used  in 
the  Reference  Case.  As  shown  in  Figure  A-9,  the  0.001  threshold  regions  of 
ADRg  in  the  two  shocks  were  reduced  slightly  for  Case  lb,  but  the  peak  ADR 
-'alues  in  the  shocks  (near  the  0.01  threshold)  increased  slightly.  The 
extent  of  the  region  of  large  disagreement  in  the  total  pressure  profile  was 
reduced  as  shown  in  Figure  A-10,  but  the  largest  disagreement  has  increased 
to  approximately  3%.  Some  noise  is  starting  to  contaminate  the  solution  near 
the  top  wall  which  is  indicated  by  the  high  ADRg  region  and  the  noise  in  the 
total  pressure  profile.  This  noise  is  probably  caused  by  the  large  cell 
aspect  ratios  in  that  region. 

Since  grid  refinement  in  the  y-direction  did  not  eliminate  the  regions  of 
high  ADR  and  total  pressure  disagreement,  an  attempt  was  made  to  refine  the 
grid  in  the  streamwise  (x)  direction.  The  upstream  (left)  boundary  was  moved 
closer  to  the  bow  shock  by  reducing  the  mesh  stretch  factor  to  the  left  of 
the  cowl  lip,  thereby  increasing  the  number  of  points  in  the  bow  shock. 
Twelve  equal  spaced  mesh  points  were  also  added  just  downstream  (to  the 
right)  of  the  bleed  slot  lip  to  increase  the  resolution  at  the  downstream  end 
of  the  slot  shock.  Figure  A-ll  shows  the  ADRg  contours  for  this  case  (Case 
2)  with  the  overlayed  grid.  This  grid  includes  refinement  of  the  grid  in  the 
streamwise  direction  and  the  same  y-direction  grid  as  was  used  for  Case  la. 
As  shown  by  the  ADRg  contours,  significant  reductions  were  made  in  the  ADRg 
greater  than  0.001  threshold  regions  downstream  of  the  slot  shock.  However, 
a  comparison  of  the  total  pressure  profiles  with  and  without  the  x-direction 
refinements  (Figure  A-12)  shows  that  the  agreement  with  the  measured  data  is 
not  improved  significantly.  Significant  reductions  in  the  noise  near  the 
upper  wall  did  occur  (indicated  by  the  reduced  ADRg  in  this  region)  since  the 
cell  aspect  ratio  was  decreased  by  the  additional  x-direction  mesh. 
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The  Case  la  results  shown  in  Figures  A-6  and  A-7  were  reviewed  again 
since  the  refinements  to  the  grid  in  the  streamwise  direction  did  not  produce 
significant  improvements  in  the  total  pressure  profile.  These  ADR  contour 
plots  indicated  a  large  region  of  high  ADR  in  the  bleed  chamber.  Therefore  a 
run  was  made  with  the  grid  in  the  bleed  chamber  increased  from  42  x  18  cells 
to  51  x  26  cells.  Color  contour  plots  from  this  run  (Case  3)  showed  a 
significant  reduction  in  ADR  in  the  bleed  chamber,  but  almost  no  change  in 
ADR  contours  in  the  inlet  region.  The  total  pressure  profile  from  this  run 
(shown  in  Figure  A-13  compared  with  Case  la)  showed  some  changes  in  the 
boundary  layer  near  the  lower  wall,  but  the  disagreement  between  y  =  0.0075 
and  0.03  meters  did  not  change. 

To  investigate  the  convergence  effects  on  ADR  contours  and  total  pressure 
profiles,  the  Reference  Case  was  restarted  and  run  an  additional  20,000 
cycles.  This  case  (Case  4)  had  then  completed  the  same  number  of  cycles  as 
was  completed  for  Case  la.  A  comparison  of  the  ADR  contours  from  the 
Reference  Case  and  Case  4  showed  that  ADR  did  not  change  significantly  with 
the  additional  cycles.  However,  as  shown  in  Figure  A-14,  the  additional 
20,000  cycles  has  resulted  in  a  more  fully  developed  boundary  layer  profile 
downstream  of  the  bleed  slot  thus  reducing  the  maximum  disagreement  in  the 
total  pressure  profile  to  about  3Z.  Figure  A-15,  which  compares  the  results 
from  Case  la  with  Case  4,  shows  that  the  effect  of  the  additional  24  y- 
direction  cells  on  the  total  pressure  profile  at  the  same  number  of  cycles  is 
to  reduce  the  maximum  disagreement  by  about  IX.  This  convergence  study  shows 
that  of  the  2%  improvement  in  total  pressure  claimed  in  Figure  A-8,  about  IX 
is  due  to  the  increased  y-direction  grid  and  the  other  IX  is  due  to 
convergence  in  this  part  of  the  flow  field. 

Examination  of  the  ADRg  contours  relative  to  mesh  density  indicated  that 
ADRg  was  high  in  regions  where  the  grid  was  too  coarse  to  properly  resolve 
the  gradients  in  the  flowfield.  Refinements  to  the  grid  in  the  direction 
normal  to  the  gradients  produced  the  greatest  reductions  in  ADRg.  In  regions 
of  the  flow  field  where  the  steep  gradients  are  diagonal  to  the  mesh, 
refinement  was  required  to  achieve  cell  aspect  ratios  close  to  unity.  If  the 
cell  aspect  ratio  was  near  unity  and  ADRg  remained  high,  refinements  to  the 
grid  were  required  in  both  directions. 

A  study  of  the  change  in  the  predicted  total  pressure  profiles  as 
compared  to  the  measured  total  pressure  profiles  showed  that  satisfying  the 
recommended  ADR  thresholds  was  not  sufficient  to  guarantee  an  accurate 
solution.  This  is  consistent  with  one  conclusion  of  Appendix  C  which  says 
that  there  is  not  a  direct  correlation  between  peak  ADR  level  and  solution 
accuracy.  The  disagreements  that  remain  could  also  be  caused  by 
uncertainties  in  the  experimental  measurements,  uncertainties  in  the 
numerical  boundary  conditions,  or  by  improper  geometry  resolution  at  both 
ends  of  the  bleed  slot. 

A  comparison  of  the  simultaneous  changes  in  ADRg  and  ADRq  for  these  grid 
changes  was  made.  In  regions  of  the  flow  field  where  the  steep  gradients 
were  diagonal  to  the  mesh  and  the  cell  aspect  ratio  was  large,  a  high  ADRg  or 
ADRq  would  indicate  that  mesh  refinement  was  needed  in  the  direction  normal 


-  A-4  - 


to  the  velocity  component  used  to  calculate  ADR.  For  example,  in  the 
expansion  region  in  Figure  A-3,  ADRp  (streamvise  momentum)  was  high 
indicating  that  normal  (y)  direction  mesh  was  needed  to  make  the  cell  aspect 
ratio  closer  to  unity.  This  effect  was  also  seen  for  high  ADRq  regions  and 
streamwise  direction  refinements.  In  regions  of  the  flow  field  where  the 
gradients  were  aligned  to  the  mesh,  all  three  ADRs  tended  to  be  near  the  same 
level,  so  direction  for  grid  refinement  could  not  be  inferred  solely  from 
ADRp  and  ADRq. 


4.  Conclusions  and  Recommendations 

The  use  of  ADR  to  guide  the  localized  grid  refinement  for  a  very  complex 
flow  analysis  has  been  shown  in  this  study.  The  regions  of  the  flow  field 
where  the  grid  was  too  coarse  to  properly  resolve  the  flow  gradients  produced 
high  values  of  ADR.  The  high  ADR  values  were  reduced  and  the  solution 
accuracy  improved  when  more  mesh  was  added. 

The  use  of  ADR  thresholds  as  a  quantitative  measure  of  solution  accuracy 
has  not  been  shown  in  this  study.  As  the  ADR  levels  were  reduced  to  the 
recommended  levels,  the  predicted  total  pressure  profiles  were  brought  into 
better  agreement  with  the  experimental  results,  however,  some  disagreement 
remained  between  the  computed  and  measured  total  pressure  profiles.  This 
conclusion  is  consistent  with  the  conclusion  of  the  Appendix  C  study  of  ADR 
for  a  ID  model  equation. 

For  this  study,  the  experience  of  the  CFD  user  was  utilized  to  determine 
that  the  high  ADRp  regions  could  be  reduced  by  adding  mesh  in  the  direction 
normal  to  the  gradient  or  in  a  way  to  reduce  the  cell  aspect  ratio.  It  would 
be  advantageous  to  develop  guidelines  using  ADR  that  could  be  given  to  novice 
CFD  users  or  could  be  programmed  into  an  adaptive  mesh  generator.  The 
investigation  into  using  ADRp  and  ADRq  to  guide  the  grid  refinement  direction 
showed  promise,  but  additional  research  is  required. 

Several  recommendations  for  continuing  work  on  this  problem  are: 

1)  Add  more  x-direction  mesh  in  the  bleed-slot  lip  shock  region. 

This  can  be  accomplished  with  the  current  mesh  generator  but 
has  not  been  performed  at  present  due  to  time  constraints. 

2)  Use  an  elliptic  mesh  generator,  preferably  with  multi-block 
capability,  to  better  define  the  geometry  at  both  ends  of  the 
bleed  slot  and  at  the  cowl  lip. 

3)  Demonstrate  the  use  of  ADR  on  other  test  problems  and  with 
other  codes. 
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Figure  A-2.  Reference  Case  -  66  Y-Diredbn  Cells 


Figure  A-3.  Reference  Case  -  66  Y-Dkeciion  Celts 
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Figure  A-5.  Total  Pressure  Comparison  -  Reference  Case  vs.  Experiment 
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Figure  AS.  Total  Pressure  Comparison  -  Reference  Case  and  Case  la. 
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Figure  A-10.  Total  Pressure  Comparison  -  Case  la.  and  1b. 
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Figure  A-12.  Total  Pressure  Comparison  -  Case  la.  and  2. 
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Figure  A-13.  Total  Pressure  Comparison  -  Case  la.  and  3. 
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Figure  A-14.  Total  Pressure  Comparison  -  Reference  Case  and  Case  4. 
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APPENDIX  B: 

Modified  Equation  Analyses  for  Spatial  Truncation  Error  Effects 
in  Numerical  Solutions  of  Burgers  Equation 

by 

GC  Paynter 


1 .  Summary 

The  contract  effort  resulted  in  the  identification  of  an  error  monitor, 
the  artificial  diffusion  ratio,  useful  for  guiding  grid  adjustment  and 
establishing  solution  accuracy.  A  follow-on  IRAD  study  was  completed  to 
explore  the  use  of  ADR  for  algorithms  other  than  MacCormack' s  explicit  and 
to  establish  a  relationship  between  the  ADR  level  and  solution  accuracy. 

The  follow-on  study  used  Burgers  equation  as  an  analog  to  the  NS  equations. 
As  part  of  this,  modified  equations  were  developed  for  the  linear  viscous 
Burgers  equation  for  both  the  MacCormack  explicit  and  the  Beam-Warming 
implicit  algorithms.  The  modified  equation  is  the  finite  difference 
equation  written  as  a  sum  of  the  partial  differential  equation  and  the 
truncation  error  with  the  lead  time  derivative  terms  in  the  truncation 
error  replaced  by  spatial  derivatives  through  a  manipulation  of  this 
equation.  The  modified  equation  portion  of  the  follow-on  study  is  reported 
below.  It  was  found  that  the  lead  truncation  error  term  with  the 
MacCormack  algorithm  is  dispersive  and  that  the  lead  term  with  the  Beam- 
Warming  algorithm  is  dissipative. 


2.  Introduction 

In  Computational  Fluid  Dynamics  (CFD),  algebraic  approximations  to  the 
Navier-Stokes  equations  (or  a  subset)  are  solved  for  the  values  of  the 
dependent  variables  at  each  point  on  a  grid  of  points  distributed  over  a 
flow  domain  of  interest.  Defining, 

PDE  =  Partial  Differential  Equation 

FDE  =  Finite  Difference  Equation 

TE  =  Truncation  Error 

Then  PDE  =  FDE  +  TE 

Since  it  is  usually  assumed  that 


FDE  =  0, 


This  implies  that  we  are  actually  solving, 

PDE  =  TE 

Solving  the  FDE  =  0  is  thus  equivalent  to  solving  the  PDE  with  a 
nonphysical  source  term  equal  to  the  truncation  error. 

Truncation  error  effects  can  be  classified  as  dissipative,  dispersive, 
or  diffusive,  as  discussed  in  Anderson,  et  al.  (1984).  The  truncation 


error  can  be  derived  from  the  FDE  by  substituting  Taylor  series  expansions 
about  a  point  for  the  dependent  variable  values  at  neighboring  points 
appearing  in  the  FDE.  Dissipation,  associated  with  the  even  derivative 
terms  in  the  truncation  error  expression,  acts  like  an  "artificial 
viscosity"  to  reduce  all  gradients  in  dependent  variables  in  the  solution. 
Dispersion,  associated  with  the  odd  derivative  terms  in  the  truncation 
error,  acts  to  distort  the  predicted  values  of  dependent  variables  in  the 
vicinity  of  wave  fronts  (like  shocks).  The  combined  effects  of  dissipation 
and  dispersion  are  referred  to  as  diffusion. 

In  addition,  it  is  common  to  introduce  explicit  "artificial  viscosity" 
into  the  solution  to  control  dispersive  error  effects.  The  problem  with 
this  is  that  adding  artificial  dissipation  to  eliminate  dispersion 
introduces  a  dissipative  error  of  unknown  magnitude  and  thus  it  has  an 
unknown  effect  on  solution  accuracy.  The  work  reported  in  this  appendix 
addresses  the  relationship  between  ADR  and  solution  accuracy  through 
modified  equation  analyses  of  the  MacCormack  and  Beam-Warming 
approximations  to  the  linear  Burgers  equation. 

The  remainder  of  this  appendix  consists  of  sections  on  Burgers 
equation,  the  modified  equation  analyses,  a  discussion  of  results,  and 
conclusions.  Numerical  solutions  to  Burgers  equation  using  the  two 
algorithms  and  comparisons  between  numerical  and  analytic  results  to 
address  ADR  as  a  measure  of  solution  accuracy,  are  reported  in  Appendix  C. 


3.  Burgers  Equation 

As  noted  in  Anderson,  et  al.  (1984)  Burgers  equation  is  a  simple 
nonlinear  analog  of  the  NS  equations.  Both  have  an  unsteady  term,  a 
convective  term  and  a  viscous  term.  The  most  important  difference  is  that 
Burgers  equation  is  a  scalar  equation  and  the  NS  equations  are  a  set  of 
vector  equations  of  the  same  form.  Burgers  equation  is  useful  for 
evalution  of  error  monitors  because  it  is  recognized  as  an  analog  to  the  NS 
equations,  analytic  solutions  are  available,  programming  for  numerical 
solutions  using  a  variety  of  algorithms  is  simple  enough  to  be  feasible, 
and  development  of  the  modified  equation  for  a  given  algorithm  is  feasible. 
The  viscous  Burgers  equation  can  be  written, 

ut+  Ex=  0  (B-l) 

where, 

2 

E  =  F  -  uu  ,  F  =  cu  +  bu 
x  2 

If  b  =  0,  the  equation  is  linear.  If  y  =  0,  the  equation  is  inviscid.  If 
b  =  0  and  y  =  0,  the  equation  becomes  the  familiar  wave  equation. 


4.  Modified  Equation  Analyses 

If  the  PDE  is  a  function  of  both  space  and  time,  the  truncation  error 
associated  with  an  FDE  approximation  of  the  PDE  is  also  a  function  of  both 
space  and  time.  The  FDE  is  an  algebraic  expression  for  the  value  of  the 
dependent  variable(s)  at  a  point  in  the  grid  at  a  new  time  level  in  terms 
of  the  value  for  the  point  at  the  old  time  level  and  values  at  neighboring 
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points  at  either  old  or  new  time  levels.  The  points  involved  in  the  FDE 
are  referred  to  as  the  finite  difference  stencil.  An  expression  for  the 
truncation  error  of  the  FDE  is  obtained  by  substituting  Taylor  series 
expansions  for  the  dependent  variable(s)  about  one  of  the  points  in  the 
stencil  for  the  values  at  the  other  points  into  the  FDE.  As  noted  in  the 
introduction,  this  results  in  an  equation  of  form, 

PDE  -  TE  *  0 

When  this  equation  is  manipulated  to  eliminate  the  lowest  order  time 
derivative  terms  appearing  in  the  truncation  error,  the  resulting  equation 
is  called  the  modified  equation.  The  modified  equation  suggests  the 
spatial  effect  of  the  lead  or  dominant  truncation  error  terms  on  the 
spatial  solution. 

Because  of  the  complexity,  it  was  impractical  to  develop  modified 
equations  for  the  nonlinear  form  of  the  viscous  Burgers  equation  in  the 
time  available  for  the  study.  Modified  equations  are  thus  reported  for 
the  linear  form  of  the  viscous  Burger's  equation,  b  =  0  in  Equation  B-l. 


4.1  The  Modified  Equation,  MacCormack  Algorithm 

The  MacCormack  explicit  algorithm,  described  in  Anderson,  et  al. 
(1984),  uses  forward-time  forward-space  differencing  for  a  predictor  step 
and  forward-time  backward-space  differencing  for  a  corrector  step.  The 
values  of  the  dependent  variables  resulting  from  the  predictor  and 
corrector  steps  are  treated  as  temporary  values  at  time  n  and  the  value  of 
the  dependent  variable  at  n+1  is  set  equal  to  the  average  of  the  values 
resulting  from  the  predictor  and  corrector  steps.  It  should  be  noted  that 
the  spatial  derivatives  in  the  corrector  step  are  computed  using  the 
temporary  values  of  the  dependent  variable  from  the  predictor  step. 

With  reference  to  Equation  B-l,  the  predictor  step  is, 


where  u*  is  a  temporary  value  of  "u"  at  n.  The  corrector  step  is, 
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where  u**  is  the  temporary  value  of  "uM  at  n  resulting  from  the  corrector 
step.  Averaging  the  temporary  values  of  "u”  resulting  from  the  predictor 
and  corrector  steps,  the  FDE  for  the  value  of  "u"  at  n+1  is  obtained. 
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Note  that  E  j  and  E  are  computed  using  the  temporary  values  of  "u" 
resulting  from  the  predictor  step. 

Selecting  point  n,j  about  which  to  expand,  Taylor  series  expansions  are 
needed  for  un+*j,  Enj+i,  E* j ,  and  E*j_j.  When  these  expansions  are 
substituted  into  the  FDE,  Equation  B-4,  an  equation  of  form  PDE  -  TE  =  0  is 
obtained.  As  reported  by  Baldwin,  et  al.  (1977),  when  this  equation  is 
manipulated  to  eliminate  the  time  derivative  terms,  the  modified  equation 
results. 

u  +  cu  -  uu  +  -§-  f  Ax2  -  c2At2  )  u  +  ...  HOT  =  0  (B-5) 

t  x  xx  6  (  J  xxx  v 

where  HOT  =  higher  order  terms 


4.2  The  Modified  Equation,  Beam-Warming  Algorithm 

From  Anderson,  et  al.  (1984),  the  Beam-Warming  algorithm  for  Equation 
B-l  can  be  written, 
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The  finite  difference  stencil  thus  involves  points  j-1,  j,  and 
j+1  at  times  n  and  n+1.  Selecting  point  n,j  about  which  to 
expand,  Taylor  series  expansions  are  needed  for  the  remaining 
points  in  the  stencil. 
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Restricting  the  analysis  to  the  linear  form  of  Equation  B-l,  the  FDE  can  be 
written, 
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Substituting  the  Taylor  expansions  from  Equation  B-7  into  Equation  B-6,  an 
equation  of  form  PDE  -  TE  =  0  is  obtained.  This  may  be  written, 
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Note  that  the  TE  =  0  (At,  Ax^). 

The  modified  equation  is  obtained  for  this  algorithm  by  manipulating 
Equation  B-8.  This  manipulation  involves  differentiating  Equation  B-8  with 
respect  to  x  and/or  t  and  multiplying  it  by  a  coefficient  such  that  when  it 
is  added  to  the  original  equation,  a  selected  time  derivative  term  is 
eliminated.  Vhen  this  is  done  in  a  sequential  fashion,  time  derivative 
terms  up  to  a  selected  order  are  replaced  in  Equation  B-8  by  equivalent 
space  derivative  terms.  Following  Anderson  (1984)  this  is  done  by 
constructing  a  table  in  which  the  equation  across  the  top  of  the  table  is 
Equation  B-8  and  the  terms  listed  on  the  left  side  are  the  manipulations 
performed  on  a  given  line  of  the  table  to  eliminate  a  time  derivative  term. 
When  the  equations  on  each  line  of  the  table  are  summed  together,  the 
modified  equation  is  the  result.  Table  B-l  shows  this  manipulation  for  the 
Beam-Warming  algorithm  applied  to  the  linear  viscous  Burgers  equation. 


5.  DISCUSSION 

With  reterence  to  the  modified  equation  for  the  MacCormack  algorithm, 
the  lead  truncation  error  term  (the  term  with  the  lowest  order  spatial 
derivative)  involves  a  third  order  spatial  derivative.  This  suggests  that 
the  dominant  truncation  error  effect  should  be  dispersive.  Implicit  in  this 
conclusion  however  is  the  difficult  to  prove  assumption  that  the  lowest 
order  term  is  the  dominant  term.  It  should  be  possible  to  test  this 
assumption  for  an  algorithm  by  constructing  analytic  solutions  to  a  model 
equation  with  a  dispersive  error  like  effect  (ringing)  near  a  steep  gradient 
region.  This  should  result  in  an  equation  of  form, 

PDE  =  SOURCE  TERM 

as  suggested  by  McDonough  (1988).  Since  the  constructed  solution  represents 
a  solution  to  the  equation, 


FDE  =  0 


This  implies  that, 

PDE  =  TE  and  that  TE  =  SOURCE  TERM 

It  should  be  possible  to  directly  evalute  the  derivatives  of  the  TE  from 
this  equation. 
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Examination  of  the  modified  equation  for  the  linear  Burgers  equation, 
approximated  using  the  MacCormack  algorithm,  reveals  that  the  lowest  order 
term  is  zero  when  the  Courant  number  is  1.  The  numerical  results  by  Mayer, 
Appendix  C,  (recall  that  Mayer  solved  the  nonlinear  Burgers  equation) 
indicate  that  the  exact  solution  is  recovered  when  the  Courant  number  is 
near  1.  The  algorithm  is  unstable  when  the  Courant  number  is  greater  than 
1. 


With  reference  to  the  modified  equation  for  the  Beam-Warming 
algorithm,  the  lead  truncation  error  term  involves  a  second  order  spatial 
derivative.  This  suggests,  somewhat  suprisingly,  that  the  dominant 
truncation  error  effect  is  dissipative  and  that  the  magnitude  of  this  term 
is  a  function  of  the  time  step.  It  should  be  noted,  however,  that  a  third 
order  term  does  exist.  Numerical  results  by  Mayer,  Appendix  C,  show  that 
when  the  grid  is  coarse  an  oscillation  of  the  solution  (a  dispersive  error 
effect)  occurs  near  a  steep  spatial  gradient.  This  suggests  that  when  the 
grid  is  coarse,  the  third  order  term  becomes  the  dominant  term  near  a  strong 
spatial  gradient. 

Artificial  dissipation  as  implemented  by  Forester  (1985),  and  Mayer, 
appendix  C,  is  a  viscous  like  term  added  to  the  FDE  to  control  a  dispersive 
error  effect  in  the  solution.  This  term  is  typically  the  product  of  an 
artifical  viscosity  and  a  second  derivative  of  a  dependent  variable.  When 
artificial  dissipation  is  used,  instead  of  solving, 

FDE  =  0 

the  equation  being  solved  is, 

FDE  =  AD 

where  AD  «  artificial  dissipation.  This  implies  that  , 

PDE  -  TE  -  AD  =  0 

Thus  AD  is  in  effect  a  source  term  added  to  the  PDE.  The  artificial 
viscosity  as  implemented  by  Forester  and  Mayer  is  proportional  to  the  second 
derivative  of  a  dependent  variable.  Thus,  unless  a  dispersive  error  is 
present  in  the  solution,  the  artificial  viscosity  and  hence  the  AD  are  zero. 
Although  AD  is  added  to  the  FDE  to  control  dispersive  error,  some  dispersive 
error  must  be  present  in  the  solution  to  activate  the  artificial  viscosity 
model. 


The  sum  of  the  TE  and  the  AD  act  as  a  source  term  in  the  PDE  being 
solved.  AD  is  at  least  a  partial  measure  of  solution  accuracy  because  in 
regions  of  high  dispersive  error  the  source  term  is  "large".  Thus  a 
solution  with  a  locally  high  AD  would  have  a  locally  high  source  term.  A 
solution  with  a  high  area  normalized  value  of  AD  would  be  expected  to  be 
less  accurate  than  one  with  a  low  value  of  this  parameter.  Unfortunately, 
other  error  mechanisms  may  be  present.  If  the  algorithm  has  intrinsic 
dissipation,  the  source  term  of  the  PDE  could  be  high  while  the  AD  is  low. 

A  second  problem  is  that  the  relationship  between  the  magnitude  of  the 
source  term  and  the  solution  accuracy  is  highly  nonlinear,  Oran  and  Boris, 
1987.  This  relationship  may  also  be  very  problem  dependent. 
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AD  would  be  expected  to  be  most  useful  for  algorithms  where  the  lead  TE 
terms  are  dispersive  and  less  useful  for  algorithms  with  intrinsic 
dissipation.  This  hypothesis  is  supported  by  the  numerical  results  of 
Appendix  C.  In  the  numerical  results,  the  relationship  between  the  local 
and  integrated  error  and  ADR  is  obviously  weaker  for  the  Beam-Warming 
algorithm  than  for  the  MacCormack  algorithm.  This  is  true  if  the  hypothesis 
is  true  because  of  the  intrinsic  dissipation  of  the  Beam-Warming  algorithm. 
It  does  not  appear  feasible  to  establish  a  relationship  between  solution 
accuracy  and  AD  that  is  not  problem  and  algorithm  dependent. 

It  should  perhaps  be  noted  that  while  addressing  algorithms  with 
intrinsic  dissipation  as  the  primary  means  of  controlling  dispersive  error 
was  not  a  goal  of  the  present  study,  the  Beam-Warming  algorithm  does  have 
intrinsic  dissipation.  It  is  clear  from  consideration  of  the  modified 
equation  for  the  Beam-Warming  algorithm  that  a  truncation  error  could  exist 
although  the  dispersive  error  is  zero  (and  hence  the  AD  needed  to  control 
it). 


6.  CONCLUSIONS 

Modified  equations  were  developed  for  the  linear  viscous  Burgers 
equation  for  both  the  MacCormack  explicit  and  the  Beam-Warming  implicit 
algorithms.  It  was  found  that  the  lead  truncation  error  term  with  the 
MacCormack  algorithm  is  dispersive  and  that  the  lead  term  with  the  Beam- 
Warming  algorithm  is  dissipative.  A  third  order  dispersive  error  term  was 
also  present  with  the  Beam-Warming  algorithm.  AD  is  in  effect  a  source  term 
added  to  the  PDE  being  solved  and  for  algorithms  with  high  dispersive  error 
it  is  thus  a  qualitative  measure  of  solution  accuracy.  AD  would  be  less 
f  •  useful  for  algorithms  with  high  intrinsic  dissipation.  It  is  clear  that  a 

quantitative  relationship  between  AD  and  solution  accuracy  would  be 
algorithm  and  problem  dependent. 
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Table  B-l  Manipulations  to  Obtain  the  Beaming-Warming  Modified  Equation 
for  the  Linear  Viscous  Burgers  Equation 


APPENDIX  C: 

Evaluation  of  the  Artificial  Diffusion  Ratio  through  Numerical 
Solutions  to  Burgers  Equation 


by 

DV  Mayer 
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1.  SUMMARY 

The  contract  explored  a  variety  of  truncation  error  monitors  as  guides  to 
grid  selection  and  solution  accuracy.  One  of  these  error  monitors,  the 
artificial  diffusion  ratio  (ADR)  seemed  to  provide  a  simple  means  of  guiding 
grid  refinement  and  establishing  solution  accuracy.  The  objective  of  the 
present  study  was  to  determine  if  ADR  is  applicable  to  other  algorithms  and 
determine  how  the  solution  accuracy  is  related  to  ADR. 

During  the  course  of  this  study  it  was  found  that 

o  solution  accuracy  is  not  simply  related  to  ADR  thresholds, 

o  ADR  is  useful  for  defining  grid  regions  that  require  refinement,  and 

o  ADR  is  less  useful  in  defining  regions  of  high  grid  error  for  the 
Briley-McDonald  algorithm  than  for  the  MacCormack  explicit  algorithm. 

It  is  recommended  that  the  usefulness  of  artificial  dissipation  for 
algorithms  with  intrinsic  smoothing  be  investigated  and  demonstrated.  The 
computer  program  developed  for  the  present  study  may  be  readily  adapted  for 
this  purpose.  These  follow-on  studies  should  include  a  modified  equation 
analysis  to  define  the  lead  truncation  error  terms  for  the  algorithm  and  model 
equation.  An  analytic  solution  could  be  used  to  directly  evaluate  these 
terms.  A  detailed  study  of  these  terms  relative  to  the  present  (or  alternate) 
definitions  of  ADR  should  provide  valuable  insight  and  guidance  to  using  ADR 
for  assessing  solution  accuracy. 

The  objective  of  the  work  presented  in  this  appendix  was  to  determine  if 
ADR  is  applicable  to  other  algorithms  and  how  the  solution  accuracy  is  related 
to  ADR.  This  appendix  presents  details  on  Burgers  equation  including  an 
analytic  solution,  finite  difference  formulations  of  the  two  solution 
algorithms,  ADR  definitions  for  each  algorithm,  and  results  of  numerical 
solutions  compared  with  ADR  levels  and  measures  of  solution  accuracy.  The 
final  section  discusses  the  major  conclusions  of  this  study  and  suggests 
promising  avenues  for  future  investigation.  The  modified  equation  analysis  of 
the  two  algorithms  is  reported  in  Appendix  B. 


2.  BURGERS  EQUATION 

Burgers  equation  is  a  simple  nonlinear  analog  of  the  NS  equations;  it 
contains  an  unsteady  term,  a  convective  term  and  a  viscous  term.  Burgers 
equation  is  useful  for  evaluation  of  error  monitors  because  it  is  recognized 
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as  an  analog  to  the  NS  equations,  analytic  solutions  are  available, 
programming  for  numerical  solutions  using  a  variety  of  algorithms  is 
straightforward,  and  development  of  the  modified  equation  for  a  given 
algorithm  is  feasible. 


2.1  Generalized  Burgers  Equation  Formulation 

The  generalized  Burgers  equation  may  be  written  in  a  form  consistent  with 
the  usur.l  representation  of  the  NS  equations  as  follows: 


3U  3E 

at  +  ax  " 

where, 

where 


(C-l) 


U  =  u,  E  =  F  -  y  |,  F  =  cu  +  ~  (C-2) 

u  is  the  velocity,  u  is  the  fluid  viscosity,  and  c  and  b  are  free 
parameters.  Note  that  if  b  =  0,  the  linearized  Burgers  equation  results 
and  if  c  =  0  and  b  =  1,  the  nonlinear  Burgers  equation  is  obtained.  If  u 
=  0  the  equation  is  inviscid  and  if  b  =  0  and  u  =  0  the  familiar  wave 
equation  is  recovered. 

1« 

2.2  Analytic  Solution 

If  the  free  parameters  b  and  c  are  set  to  0.5  and  -1,  respectively, 
the  generalized  Burgers  equation  has  the  stationary  solution  (Anderson  et 
al.,  1984) 


ua  =  -  H1 


c(x  -  X 


tanh 


2u 


-] 


(C— 3 ) 


A  plot  of  the  analytic  solution  is  shown  in  Figure  C-l  for  several  values 
of  y.  Note  that  this  forms  a  stationary  wave  that  approximates  a 
discontinuity  (e.g.,  a  shock  wave).  This  analytic  solution  was  therefore 
felt  to  be  a  good  choice  for  investigations  of  the  error  monitor 
identified  during  the  contract. 


3.  SOLUTION  ALGORITHMS 

The  MacCormack  explicit  and  the  Briley-McDonald  implicit  algorithms 
were  chosen  for  the  present  study.  The  MacCormack  explicit  algorithm  was 
selected  because  the  2D  NS  code  used  in  the  contract  work  is  based  on  the 
MacCormack  method.  The  Briley-McDonald  algorithm  was  selected  since  it  is 
equivalent  to  the  popular  Beam-Varming  method  when  applied  to  Burgers 


C-2 


equation.  The  numerical  results  from  solution  of  Burgers  equation  using 
these  two  algorithms,  provides  a  link  between  the  results  of  the  AFOSR 
contract  work  and  extensions  to  currently  popular  Navier-Stokes  solvers 
such  as  the  ARC  (Pulliam,  1984)  and  PARC  (Cooper,  1987)  series  of  codes. 

The  next  two  sections  discuss  these  algorithms  in  more  detail  and 
provide  the  mathematical  formulations  used  in  the  ID  Burgers  equation 
studies. 


3.1  MacCormack  Explicit  Method 

The  finite-difference  expression  for  MacCormack' s  explicit  predictor- 
corrector  method  with  a  "product"  smoother  is  given  by  (Anderson  et  al., 
1984) 


-  £  Ki  -  -  <sh  -  sM  <c-4> 

«r‘  -  H0?  *  ui  -  &  K  -  o  - <s*  -  *;.,>])  <c-5) 

where  S  represents  the  artificial  dissipation  flux  at  a  point.  For  the  NS 
equations  S  is  defined  as 


|pn  -  2P°  +  Pn  1 

s"  *  G  (  |u"|  +  a")  i-iil - i - izll  '"n 

1  1  1  (Pn  +  2Pn  +  Pn  ) 

^i+1  i  +  i-l' 


(U*i‘  -  U^) 


S.  =  e  (  lUj  I  +  ai) 


i  *  _  *  * 

P.  ,  -  2P.  +  P. 

1  l+l  l 


i-1 


★  ★  ★ 
<Pi+l  +  2Pi  +  Pi-1> 


*  * 
<Ui+l  -  Ui> 


and  0  <  e  <0.5 


(C-6) 


(C— 7 ) 


Note  that  when  solving  Burgers  equation,  the  pressure,  P,  is  not  a 
dependent  variable  and  the  speed  of  sound,  a,  is  not  a  solution  parameter 
Therefore,  for  this  study,  velocity  was  used  in  place  of  pressure  and  the 
wave  speed,  c,  was  used  in  place  of  a  in  the  definition  of  the  artificial 
dissipation  flux. 


An  additive  artificial  viscosity  is  thus  defined  as 
ua,i  ‘  E  (to)  r  K.l  -  2ui  *  Vll 
where  r  = 


(C-8) 


(  |u.  |  +  c  ) 


(u.+l  +  2u.  +  ui_1) 
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The  Ax  is  included  in  the  definition  of  the  artificial  viscosity  so  that 
S  »«  ya  3u/3x  and  is  comparable  to  the  viscous  term  u  3u/3x  contained  in 
E.  Also  note  that  the  sign  on  the  smoothing  terms  in  Equations  C-4  and  C- 
5  have  been  selected  so  that  they  have  the  same  sign  as  the  viscous  term 
in  E.  This  smoother  is  triggered  by  oscillations  in  a  numerical  solution 
(a  third  order  or  dispersive  effect).  With  an  appropriate  setting  for  e 
the  dissipation  added  by  the  smoother  will  damp  the  third  order  dispersion 
error. 

For  this  study  u  is  <  0  for  certain  regions  of  the  solution  domain 
(refer  to  Figure  C-l).  This  condition  results  in  negative  or  zero  values 
for  the  artificial  viscosity  computed  from  Equation  C-8  and  causes 
numerical  difficulties.  This  problem  is  not  normally  experienced  when 
using  MacCormack's  smoother  for  the  NS  equations  since  pressure  should  be 
greater  than  zero  throughout  the  flow  field.  To  overcome  this  difficulty 
r  was  assumed  equal  to  1  for  all  conditions.  An  order  of  magnitude 
analysis  of  r  indicates  this  assumption  is  approximately  correct  if  the 
solution  is  shifted  far  from  the  origin.  For  example,  if  r  is  evaluated 
in  a  region  of  the  flow  where  u  =  100  then  r  *  0.25. 


3.2  Briley-McDonald  Implicit  Method 

The  Brilay-McDonald  method  is  an  implicit  algorithm  with  the  following 
finite-difference  expression  for  the  generalized  Burgers  equation 
(Anderson  et  al.,  1984) 

n+1 
u. 

_i _ _ 

At 


where 

A  =  ^  =  c  +  bu  (C-10) 

This  method  is  equivalent  to  the  Beam-Warming  method  when  applied  to  the 
Burgers  equation.  Smoothing  was  added  to  this  algorithm  by  replacing  u 
with  u  +  ua  where  ua  is  defined  by  Equation  C-8.  To  facilitate  the 
programming  of  the  smoother  term,  lagged  values  of  u  were  used  in  the 
computation  of  ua. 
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4.  ERROR  MONITOR  DEFINITION 

The  contract  results  indicated  that  the  magnitude  of  the  smoothing  term 
for  the  MacCormack  method  may  be  indicative  of  the  error  magnitude.  This 
may  be  true  since 

o  the  lead  truncation  error  term  involves  a  third  order  spatial 
derivative  and  is  dispersive  in  nature  (Appendix  B),  and 

o  the  MacCormack  smoother  adds  sufficient  artificial  dissipation  to 
damp  the  dispersive  effect  (i.e.,  ringing  or  oscillatory  solutions). 

Therefore,  the  magnitude  of  the  smoothing  term  may  be  of  the  same  order  of 
magnitude  as  the  lead  truncation  error  term.  For  the  contract  effort  the 
error  monitor  was  defined  as  the  ratio  of  the  smoothing  flux  to  the  total 
flux  through  computational  cell  faces  (Forester  and  Tjonneland,  1988). 

This  error  measure  was  called  the  artificial  diffusion  ratio,  ADR.  The 
following  sections  provide  mathematical  details  on  computing  ADR  for  the 
two  solution  algorithms  chosen  for  this  study. 
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4.1  ADR  for  the  MacCormack  Method 


The  following  definition  of  the  artificial  diffusion  ratio,  ADR,  was 
used  during  this  study  and  is  consistent  with  the  definition  used  in  the 
AFOSR  contract  report. 


ADR 


(C-U) 


Note  that  the  magnitude  of  the  added  artificial  smoothing  term  is  compared 
with  the  3E/3x  term  such  that  ADR  must  range  between  0  and  1.  This 
definition  is  for  the  predictor  step.  A  similar  definition  may  be 
expressed  for  the  corrector  step,  but  studies  have  shown  that  both 
definitions  have  essentially  the  same  characteristics.  ADR  based  on  the 
predictor  step  was  used  for  the  AFOSR  contract  work. 


4.2  ADR  for  the  Briley-McDonald  Method 

The  definition  of  ADR  used  for  the  Briley-McDonald  method  is 
mathematically  more  complex  than  the  definition  used  for  MacCormack' s 
method.  It  is  comparable  from  the  standpoint  that  the  magnitude  of  the 
artificial  smoothing  term  is  compared  with  the  3E/3x  term.  Defining 
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This  definition  is  equivalent  to  the  definition  of  ADR  used  for  the 
MacCormack  method. 


5.  RESULTS 

The  methods  and  the  definitions  of  ADR  discussed  in  the  previous 
sections  were  programmed  into  a  Fortran  code  that  can  be  rapidly  executed 
on  a  Digital  Equipment  Corporation  VAX  computer.  Graphics  files  were 
generated  and  procedures  for  automatically  generating  plots  were  developed. 
This  level  of  automation  allowed  rapid  investigation  of  ADR  for  the  two 
methods . 

The  first  step  was  to  determine  grid  spacing  and  fluid  properties  that 
would  result  in  a  noisy  solution  (and  thereby  trigger  the  smoother).  To 
avoid  stability  issues,  the  following  stability  criteria  suggested  by 
Tannehill  (Anderson  et  al.,  1984)  was  satisfied  for  all  numerical  studies 
presented: 


At  < 


T  ,c^)2 

|A  |  Ax  +  2y 


(C-16) 


The  solution  for  three  values  of  viscosity  (y  =  1.,  0.1,  and  0.01)  was 
computed  with  the  MacCormack  method.  Time  and  space  steps  of  0.1  and  0.5, 
respectively,  were  used  for  these  calculations.  The  numerical  and  analytic 
solution  results  are  shown  in  Figure  C-l.  The  largest  value  of  viscosity 
did  not  result  in  any  noise.  The  intermediate  viscosity  value  resulted  in 
negligible  noise.  The  smallest  value  of  viscosity  triggered  appreciable 
noise. 
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The  solution  was  next  recomputed  with  the  Briley-McDonald  method  to 
check  that  the  smallest  value  of  viscosity  would  trigger  noise  with  this 
method.  This  result  is  shown  in  Figure  C-2.  Since  the  smallest  value  of 
viscosity  triggered  noise  for  both  methods  this  value  was  used  for  the 
remainder  of  the  study. 

A  simple  grid  refinement  study  was  next  performed  for  both  algorithms. 
The  cases  were  run  first  with  the  smoother  turned  off  and  then  with  the 
smoother  turned  on.  For  the  cases  with  smoothing  the  smoothing 
coefficient,  e,  was  fixed  at  a  value  of  0.25.  The  grid  spacing  was 
started  at  Ax  =  2.  and  was  reduced  by  factors  of  2  until  a  grid  spacing  of 
Ax  =  0.03125  was  reached.  Thus,  the  solution  was  computed  on  seven  grid 
levels.  For  brevity  detailed  results  will  only  be  shown  for  grid  spacings 
of  Ax  *  1.,  0.25,  and  0.0625  (i.e.,  representative  coarse,  medium,  and 
fine  grid  results).  To  avoid  stability  problems  it  was  necessary  to  reduce 
the  step  size  for  two  of  the  grids  (Ax  =  0.0625  and  0.03125).  For 
consistency  all  the  solutions  were  run  to  t  =  10. 

During  these  calculations  several  measures  of  solution  error  were 
collected.  The  objective  was  to  look  at  these  error  measures  relative  to 
ADR  and  determine  if  a  correlation  between  error  and  ADR  exists.  The 
values  collected  were  the  terror,  ADR,  Eerror,  EADR,  peak  error,  and  peak 
ADR  defined  as 


terror  =  100  - - -  (17) 

u  -  u  . 
max  min 

where  ua,  umax,  and  umjn  are  based  on  the  analytic  solution, 
x 

.  max 

Eerror  =1  |u  -  ua|  dx  (18) 

Xmin 

x 

r  max 

EADR  =  ADR  dx  (19) 

Jx  j 
min 

the  peak  error  and  peak  ADR  are  simply  the  maximum  values  for  the  absolute 
error  and  ADR  computed  for  a  given  grid  density. 

The  details  of  these  results  for  each  algorithm  are  summarized  in  the 
next  two  sections. 


5.1  MacCormack  Explicit  Results 
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The  MacCormack  explicit  results,  without  smoothing,  are  shown  in 
Figures  C-3  through  C-5.  These  figures  show  the  results  in  the  sequence  of 
coarse,  medium  and  fine  grids,  respectively.  The  top  plot  in  each  figure 
compares  the  analytic  and  numerical  solutions  for  the  velocity.  The  middle 
plot  shows  the  terror.  The  bottom  plot  shows  the  ADR  values  which  are 
identically  zero  when  smoothing  is  not  used.  As  expected  the  range  and 
extent  of  the  error  is  minimized  with  improved  grid  density. 

The  results  with  the  smoother  turned  on  are  shown  in  Figures  C-6 
through  C-8  (in  the  same  sequence  as  Figures  C-3  through  C-5).  The  extent 
of  the  high  error  region  is  reduced  with  improved  grid  density  and  the  ADR 
values  exhibit  a  similar  behavior.  The  range  of  the  error  seems  to  be 
nearly  the  same  for  the  coarse  and  medium  grids,  and  increases  for  the  fine 
grid.  The  range  of  the  ADR  values  increases  with  grid  refinement. 

Figure  C-9  shows  the  peak  error  (with  and  without  smoothing)  and  the 
peak  ADR  (with  smoothing  only)  as  a  function  of  grid  spacing  for  all  the 
grids  investigated.  Figure  C-10  shows  the  Eerror  and  EADR  as  a  function 
of  grid  spacing.  Two  features  may  be  noted  from  these  results. 

1)  There  is  a  relationship  between  the  peak  error  and  peak  ADR,  but  it 
is  not  a  simple  one-to-one  correlation  (compare  the  middle  and 
bottom  plots  in  Figure  C-9).  This  means  the  threshold  idea  proposed 
during  the  contract  is  not  supported  by  these  results.  That  is,  the 
magnitude  of  ADR  is  not  directly  related  to  the  magnitude  of  the 
error  in  the  numerical  solution. 

2)  It  is  apparent  that  the  use  of  smoothing  is  beneficial  (compare  the 
top  and  middle  plots  in  Figures  C-9  and  C-10).  The  accuracy 
(whether  measured  as  peak  error  or  Eerror)  is  significantly 
improved  for  coarse  grids.  Conversely,  the  use  of  the  smoother  does 
not  significantly  impair  the  result  for  fine  grids.  The  low  error 
observed  for  the  coarsest  grid  (Ax  =  2),  without  smoothing,  is 
probably  due  to  a  Courant  number  and  diffusion  number  effect.  For 
these  conditions  the  truncation  error  is  apparently  negligible  and  a 
nearly  exact  solution  is  recovered.  This  should  be  checked  by 
performing  a  modified  equation  analysis  for  the  nonlinear  Burgers 
equation. 


5.2  Briley-McDonald  Implicit  Results 

The  Briley-McDonald  implicit  results,  without  smoothing,  are  shown  in 
Figures  C-ll  through  C-13.  These  figures  show  the  results  in  the  sequence 
of  coarse,  medium  and  fine  grids,  respectively.  The  top  plot  in  each 
figure  compares  the  analytic  and  numerical  solutions  for  the  velocity.  The 
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middle  plot  shows  the  terror.  The  bottom  plot  shows  the  ADR  values  which 
are  identically  zero  when  smoothing  is  not  used.  As  expected  the  range  and 
extent  of  the  error  is  minimized  with  improved  grid  density.  This  trend  is 
consistent  with  the  trend  noted  for  the  MacCormack  algorithm. 

The  results  with  the  smoother  turned  on  are  shown  in  Figures  C-14 
through  C-16  (in  the  same  sequence  as  Figures  C-ll  through  C-13).  The 
extent  of  the  high  error  region  is  reduced  with  improved  grid  density  and 
the  ADR  values  exhibit  a  similar  behavior.  The  results  differ  somewhat 
from  the  MacCormack  results  since  the  extent  of  the  high  ADR  regions  are 
greater  for  the  Briley-McDonald  method.  The  range  of  the  error  increases 
with  grid  refinement  (differing  from  the  trend  noted  for  the  MacCormack 
method  where  the  coarse  and  medium  grid  results  have  nearly  the  same  range 
of  error).  The  range  of  the  ADR  values  is  nearly  the  same  for  the  coarse 
and  medium  grids,  and  increases  for  the  fine  grid  (differing  from  the  trend 
noted  for  the  MacCormack  method  where  the  range  of  ADR  values  increased 
with  grid  refinement). 

Figure  C-17  shows  the  peak  error  (with  and  without  smoothing)  and  the 
peak  ADR  (with  smoothing  only)  as  a  function  of  grid  spacing  for  all  the 
grids  investigated.  Figure  C-18  shows  the  Terror  and  EADR  as  a  function 
of  grid  spacing.  The  following  features  may  be  noted  from  these  results: 

1)  A  simple  relationship  between  the  peak  error  and  peak  ADR  does  not 
exist  (compare  the  middle  and  bottom  plots  in  Figure  C-17).  That 
is,  the  magnitude  of  ADR  is  not  directly  related  to  the  magnitude  of 
the  error  in  the  numerical  solution.  Therefore,  an  ADR  threshold  is 
not  supported  for  either  the  MacCormack  or  Briley-McDonald  method. 

2)  A  weaker  relationship  exists  between  the  peak  ADR  and  the  peak  error 
for  the  Briley-McDonald  algorithm  than  is  observed  for  the 
MacCormack  algorithm  (compare  the  middle  and  bottom  plots  of  Figures 
C-9  and  C-17).  Similarly,  the  relationship  is  weaker  between  the 
EADR  and  Terror  (compare  the  middle  and  bottom  plots  of  Figures  C- 
10  and  C-18). 

3)  It  is  apparent  that  the  use  of  smoothing  is  beneficial  (compare  the 
top  and  middle  plots  in  Figure  C-17  and  Figure  C-18).  The  accuracy 
(whether  measured  as  peak  error  or  Terror)  is  significantly 
improved  for  coarse  grids.  Conversely,  the  use  of  the  smoother  does 
not  significantly  impair  the  result  for  fine  grids.  The  low  error 
observed  for  the  coarsest  grid  (Ax  =  2),  without  smoothing,  is 
probably  due  to  a  Courant  number  and  diffusion  number  effect.  For 
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these  conditions  the  truncation  error  is  apparently  negligible  and  a 
nearly  exact  solution  is  recovered.  This  should  be  checked  by 
performing  a  modified  equation  analysis  for  the  nonlinear  Burgers 
equation. 


6.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  summary,  the  following  conclusions  have  been  drawn  from  the  1-D 
Burgers  equation  results: 

o  The  ADR  threshold  as  a  measure  of  solution  accuracy  has  not  been 
established.  There  does  exist  a  correlation  between  ADR  and  solution 
accuracy,  but  the  correlation  is  more  complex  than  a  simple  one-to- 
one  relationship.  Achieving  a  given  level  of  ADR  does  not  guarantee 
a  given  accuracy  level. 

o  The  use  of  ADR  to  pinpoint  grid  problems  has  been  demonstrated. 
Regions  of  the  flow  where  the  grid  is  too  coarse  to  accurately 
capture  gradients  results  in  noise  or  ringing  in  the  solution.  This 
noise  is  dispersive  in  nature,  activates  MacCormack  type  smoothers, 
and  results  in  correspondingly  high  values  of  ADR.  Therefore, 
regions  of  the  flow  that  have  relatively  high  levels  of  ADR  are 
candidates  for  further  grid  refinement. 

o  ADR  does  not  pinpoint  regions  of  grid  error  as  well  for  the  Briley- 
McDonald  method.  The  regions  of  high  ADR  cover  a  greater  extent  than 
the  regions  of  relatively  high  solution  error.  This  means  that  ADR 
for  the  Briley-McDonald  method  exaggerates  the  extent  of  regions  that 
require  grid  refinement. 

This  investigation  could  be  readily  extended  to  other  algorithms  using 
techniques  similar  to  those  describe  in  this  coordination  sheet. 

The  usefulness  of  ADR  for  solution  algorithms  that  use  intrinsic 
smoothing  (i.e.,  not  added  explicitly)  has  not  been  demonstrated.  The 
computer  program  written  for  the  present  study  may  be  readily  adapted  to 
explore  the  use  of  ADR  for  these  and  other  algorithms. 

The  results  of  the  modified  equation  analysis  of  the  MacCormack 
explicit  and  Briley-McDonald  implicit  algorithm  (Appendix  B)  makes  it 
feasible  to  directly  evaluate  the  lead  truncation  error  terms  for  an 
analytic  solution.  A  detailed  study  of  these  terms  relative  to  the  present 
(or  alternate)  definitions  of  ADR  may  provide  valuable  insight  and  guidance 
to  using  ADR  for  assessing  accuracy. 
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MocCormock  explicit 
Ot  -  0.100,  DX  -  0  500 
Mu  *  1.000,  Epsilon  =  0.00 
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MocCormock  explicit 
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MocCormock  explicit 
Dt  -  0.100,  DX  -  0.500 
Mu  “  0.010,  Epsilon  m  0.00 
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Figure  C-1.  MacCormack  Explicit  Results  Compared  With  Analytic  Velocity  Profiles 
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MocCormock  explicit 
Dt  -  O.IOO,  DX  -  0  250 
Mu  ■  0.010,  Epsilon  *»  0.00 
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MocCormock  explicit 
Dt  ■  0.100,  DX  -  0.250 
Mu  ~  0.010,  Epsilon  m  0.00 
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Figure  C-4.  MacCormack  Explicit  Results:  Medium  Grid,  No  Smoothing 
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MacCormock  explicit 
Dt  -  O.IOO,  DX  -  1.000 
Mu  =  0.010,  Epsilori  =  0.25 
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MacCormock  explicit 
Dt  -  0.100,  DX  -  1.000 
Mu  =»  0.010,  Epsilon  -  025 


Figure  C-6.  MacCormack  Explicit  Results:  Coarse  Grid,  With  Smoothing 
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WocCormock  explicit 
Mu  “  0.010,  Epsilon  »  0.00 


1.0 


OX 


WocCormock  explicit 
Mu  —  0.010,  Epsilon  -  0.25 


OX 


Figure  C-10.  MacCormack  Explicit  Results:  Eerrorand  ZADR 
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Briloy— McDortold  implicit 
Dt  -  0.100,  DX  -  1.000 
Mu  =  0.010,  Epsilon  =  0.25 


Figure  C-1 4.  Briley-McDonald  Implicit  Results:  Coarse  Grid,  With  Smoothing 
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Figure  C-18.  Briley-McDonald  Implicit  Results:  Zerror  and  ZADR 


C-28 


APPENDIX  D: 

Source  Code  Listing  for  Burgers  Equation  Study 
PROGRAM  BURGER 

C**  Author:  Dave  Mayer 

C**  Boeing  Advanced  Systems 

C**  PO  Box  3703,  M/S  33-14 

C**  Seattle,  VA  98124 

C**  (206)  241-4403 

INCLUDE  'BURGER. INC/ LI ST' 

IMPLICIT  REAL*8  (A-H,0-Z) 

PARAMETER  (MAXN=500+1,  MAXI = 500+ 1 ) 

COMMON  /IOUNITS/ 

$  LUE,  LUI,  LUO,  LUP 

COMMON  /VAP. / 

$  TITLE,  I SCHEME 

$  ,  X(MAXI) ,  XMIN,  DX,  NDX,  IMAX,  ISKIP,  XO 

$  ,  T(MAXN) ,  TMIN,  DT,  NDT,  NMAX,  NSKIP 

$  ,  VIS,  EPSILON 

$  ,  ADR1 ( MAXI , MAXN ) ,  ADR2 ( MAXI , MAXN ) ,  ADR3 ( MAXI , MAXN ) 

$  ,  F(MAXI),  A (MAXI ) 

$  ,  PHI (MAXI, MAXN),  PHIB ( MAXI , MAXN) 

$  ,  B,  C,  PHIA( MAXI, MAXN) 

$  ,  PHII ,  PHIS1 ,  PHIS2 

COMMON  /STATS/ 

$  AMERR(MAXN) ,  SIGMA(MAXN),  RMSERR(MAXN) 

S  ,  SUMPHI(MAXN) ,  SUMERR(MAXN) 

$  ,  SUMADRl(MAXN) ,  SUMADR2(MAXN) ,  SUMADR3(MAXN) 

$  ,  PEAKERR(MAXN) 

$  ,  PEAKADR1 ( MAXN ) ,  PEARADR2 ( MAXN ) ,  PEARADR3 ( MAXN ) 

CHARACTERS  TITLE 


CALL  10 

100  CALL  INPUT 
CALL  INIT 
CALL  ANALYTIC 
IF( I SCHEME  .EQ.  1)  THEN 

C**  Solve  using  MacCormack  explicit  procedure 
CALL  MAC 

ELSE  IF (I SCHEME  .EQ.  2)  THEN 

C**  Solve  using  the  Briley-McDonald  implicit  procedure 
CALL  BRILEY 
ENDIF 

CALL  STATS 
CALL  EGG 
CALL  EGGS 
CALL  uP 
GO  TO  100 
END 
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SUBROUTINE  ANALYTIC 
C**  Compute  the  analytic  solution. 

INCLUDE  'BURGER. INC' 

IF(B  .EQ.  -1.  .AND.  C  .EQ.  0.5)  THEN 
C**  Compute  analytic  solution 

PI  =  AC0S(  -1.  ) 

DO  300  N*1,NMAX,NSKIP 
DO  200  1*1 , IMAX, ISKIP 
PHIA(I,N)  =  -C  /  B 

$  *  (1.  +  TANH(  C  *  (  X(I)  -  XO  )  /  (2.  *  VIS)  )  ) 

200  CONTINUE 

300  CONTINUE 

ENDIF 
RETURN 
END 
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SUBROUTINE  BRILEY 

C**  Solve  using  the  Brile^-McDonald  implicit  solution  procedure 


C** 


C** 


100 

C** 


C** 


200 

C** 


300 

C** 


400 


INCLUDE  'BURGER. INC' 

DIMENSION  AA(MAXI),  BB(MAXI) ,  CC(MAXI),  DD(MAXI) 

DO  400  N=1,NMAX-1 

Load  coefficients  for  boundary  at  1=1  (Dirichlet  condition) 
DD(1)  =  1. 

AA(1)  =  0. 

CC(1)  =  PHIS1 
Compute  A  and  F 
DO  100  1=1 , IMAX 

A(I)  -  C  +  B  *  PHI(I ,N) 

F(I)  =  C  *  PHI(I,N)  +  B  *  PHI(I ,N)**2  /  2. 

CONTINUE 

Load  coefficients  for  interior  points 
CN  =  DT  /  DX 
DO  200  1=2 , IMAX-1 

Note  that  the  smoothing  viscosity  is  based  on  lagged  values 

DN  =  (  VIS  +  VISA(I,N)  )  *  DT  /  DX**2 

BB(I)  =  -  CN  /  2.  *  A(I-l)  -  DN 

DD(I)  =  1.  +  2.  *  DN 

AA(I)  =  CN  /  2.  *  A(I+1)  -  DN 

CC(I)  =  (  -  CN  /  2.  *  A(I-l)  )  *  PHI(I-1,N) 

$  +  PHI(I,N) 

$  +  (  CN  /  2.  *  A(I+1)  )  *  PHI(I+1 ,N) 

$  -  CN  /  2.  *  (  F(I+1)  -  F(I-l)  ) 

CONTINUE 


$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


Load  coefficients  for  boundary  at  I=IMAX  (Dirichlet  condition) 
BB(IMAX)  =  0. 

DD(IMAX)  =  1. 

CC(IMAX)  =  PHIS2 

CALL  TRIDIAG ( 1 , IMAX , BB , DD , AA , CC ) 

DO  300  1=1, IMAX 
PHI(I,N+1)  =  CC(I) 

CONTINUE 

DO  400  1=2, IMAX-1 


Compute  ADR 
FLUXC  =  0.5  *  (  F(I+1) 
+  A(I+1)  * 
-  A(I-l)  * 


-  F(I-l) 

(  PHI(I+1,N+1)  -  PHI(I+1 ,N)  ) 

(  PHI(I-1,N+1)  -  PHI(I-1,N)  )  ) 


FLUXV  =  VIS 

*  (  PHI(I+1 ,N+1 )  -  2.*PHI(I,N+1)  +  PHI(I-1 ,N+1)  )  /  DX 
FLUXS  =  VISA(I,N) 

*  (  PHI(I+1,N+1)  -  2. *PHI(I ,N+1)  +  PHI(I-1,N+1)  )  /  DX 
ADR1(I,N+1)  =  ABS(  FLUXS  ) 

/  (  ABS(  FLUXC  -  FLUXV  )  +  ABS(  FLUXS  )  +  l.E-24  ) 
ADR2(I,N+1)  =  ABS(  FLUXS  ) 

/  (  ABS(  FLUXC  )  +  ABS(  FLUXV  )  +  ABS(  FLUXS  ) 

+  l.E-24  ) 

ADR3(I ,N+1)  =  ABS(  FLUXS  ) 

/  (  ABS(  FLUXC  )  +  ABS(  FLUXS  )  +  l.E-24  ) 


CONTINUE 

RETURN 

END 
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SUBROUTINE  EGG 

C**  Generate  plot  file  for  the  Boeing  Engineering  Graphics 
C**  Generator  program. 

INCLUDE  'BURGER. INC' 

WRITE (LUE, 8000) 

8000  F0RMAT('(F5.0,7E11.0)') 

WRITE ( LUE, 8005)  TITLE,  DT,  DX 
$  ,  VIS,  EPSILON 

8005  FORMAT ( ' *DUPT' / 

$  ' $' A60/ 

$  ' $Dt  ='F6.3' ,  DX  «'F6.3/ 

$  ' $Mu  =' F6. 3' ,  Epsilon  ='F6.2) 

DO  200  N=1 ,NMAX,NSKIP 
K  =  K  +  1 

WRITE (LUE, 8010)  DX,  SUMERR(N) ,  SUMADRl(N) 

$  ,  PEAKERR(N),  PEAKADRl(N) 

8010  FORMAT ( ' *PAP  DS'1PE11.3/ 

$  '*PAR  SUMERR'E11.3/ 

$  ' *PAR  SUMADR1'E11.3/ 

$  ' *PAR  PEAKERR' Ell . 3/ 

$  ' *PAR  PEAKADR1'E11.3) 

WRITE (LUE, 8020)  K 
8020  FORMAT ( 'RUN' 12) 

WRITE (LUE, 8030) 

8030  FORMAT (T5 ' I ' T9 ' X ' T20 ' Ua ' T3 1 ' U ' T42 ' *Er ror ' 

$  T53 ' ADR1 ' T64 ' ADR2 ' T7  5 ' ADR3 ' ) 

DO  100  1=1 , IMAX, ISKIP 
PHIMIN  =  PHI SI 
PHIMAX  =  PHIS2 

PCTERR  =  100.  *  (  PHI(I,N)  -  PHIA(I,N)  )  /  (  PHIMAX  -  PHIMIN  ) 
WRITE (LUE, 8040)  I,  X(I),  PHIA(I,N),  PHI(I,N) 

$  ,  PCTERR,  ADR1(I,N) ,  ADR2(I,N),  ADR3fI,N) 

8040  FORMAT(I5, 1P7E11 .3) 

100  CONTINUE 

WRITE (LUE, 8050) 

8050  FORMAT ( '*E0F' ) 

200  CONTINUE 
RETURN 
END 
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SUBROUTINE  EGGS 

C**  Generate  plot  file  of  statistical  error  measures 

C**  for  use  with  the  Boeing  Engineering  Graphics  Generator  program. 

INCLUDE  'BURGER. INC' 

LU  =  9 
K  .  K  +  1 

IF(  K  .EQ.  1  )  THEN 
WRITE ( LU , 8000 ) 

8000  FORMATS (10E11.0)') 

WRITE (LU, 8010)  TITLE,  VIS,  EPSILON 
8010  FORMAT ( '$'A60/ 

$  ' $Mu  ='F6.3' ,  Epsilon  ='F6.2) 

WRITE (LU, 8020) 

8020  FORMAT ('RUN  1') 

WRITE (LU, 8030) 

8030  FORMAT (T4'DX'T15'SUMPHI'T26' SUMERR' 

$  T37 ' SUMADR1 ' T48 ' SUMADR2 ' T59 ' SUMADR3 ' 

$  T70' PEAKERR' 

$  T81 ' PEAKADR1 ' T9 1 ' PEAKADR2 ' T102 ' PEAKADR3 ' ) 

ENDIF 

WRITE (LU, 8040)  DX,  SUMPHI(NMAX) ,  SUMERR(NMAX) 

$  ,  SUMADRl(NMAX) ,  SUMADR2(NMAX) ,  SUMADR3 (NMAX ) 

$  ,  PEAKERR (NMAX) 

$  ,  PEAKADR1 (NMAX) ,  PEAKADR2 ( NMAX ) ,  PEAKADR3 (NMAX) 

8040  F0RMAT(1P10E11 . 3) 

RETURN 

END 


SUBROUTINE  INIT 
C**  Initialize  variables 


INCLUDE  'BURGER. INC' 

C**  Set  indices 

NMAX  =  NDT  -t-  1 
IMAX  =  NDX  +  1 

C**  Set  time  and  space  grid  coordinates 
T(l)  =  TMIN 
DO  100  N=2,NMAX 
T(N)  =  T(N-l)  +  DT 
100  CONTINUE 

X(l)  =  XMIN 
DO  200  1=2, IMAX 
X(I)  =  X(I-l)  +  DX 
200  CONTINUE 

C**  Set  initial  conditions 
DO  300  1=1, IMAX 
PHI (1,1)  =  -C  /  B 

$  *  (1.  +  TANH(  C  *  (  X(I)  -  XO  )  /  (2.  *  VIS)  )  ) 

300  CONTINUE 

C**  Set  boundary  conditions 
PHI SI  =  PHI( 1 , 1 ) 

PHIS2  =  PHI ( IMAX, 1) 

RETURN 

END 
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SUBROUTINE  INPUT 
C**  Read  input  data 

INCLUDE' BURGER. INC' 

C**  Input  values  for 
C**  Title, 

C**  Ischeme:  1  — >  MacCormack  explicit 

C**  2  — >  Briley-MacDonald  implicit 

C**  Minimum  value  for  X,  Delta  X,  Number  of  Delta  X,  Output  skip  factor 

C**  Minimum  value  for  t,  Delta  t,  Number  of  Delta  t,  Output  skip  factor 

C**  Output  skip  factor 

C**  Viscosity,  Artificial  viscosity  coefficient 

C**  B,  C  — >  The  coefficients  in  the  generalized  Burger's  eq. 

C**  Ut  +  (C  +  B  *  U)  *  Ux  =  VIS  *  Uxx 

READ(LUI , 5000, END= 100)  TITLE 
5000  FORMAT (A60) 

READ(LUI,*)  ISCHEME 
READ(LUI,*)  XMIN,  DX,  NDX,  I SKIP 
XMAX  =  XMIN  +  NDX  *  DX 
XO  =  (XMAX  +  XMIN)  /  2. 

READ(LUI,*)  TMIN,  DT,  NDT,  NSKIP 
TMAX  =  TMIN  +  NDT  *  DT 
READ(LUI,*)  VIS,  EPSILON 
READ(LUI,*)  B,  C 
RETURN 

100  STOP  'End  of  input  data  detected' 

END 
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SUBROUTINE  10 

C**  Initialize  logical  units  and  files  for  Input/Output 
INCLUDE  'BURGER. INC' 

LUI  =  5 

C  0PEN(UNIT=LUI ,  FI LE=' BURGER. DAT ' ,  STATUS= ' UNKNOWN ' ) 

LUO  =  6 

C  0PEN(UNIT=LU0,  FILE* 'BURGER. OUT ' ,  STATUS* ' UNKNOWN' ) 

LUP  *  7 

C  OPENfUNIT-LUP,  FILE* ' BURGER. PLT' ,  STATUS* ' UNKNOWN ' ) 

LUE  *  8 

C  OPEN(UNIT=LUE,  FILE* ' BURGER. GGP ' ,  STATUS* ' UNKNOWN ' ) 

RETURN 
END 
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SUBROUTINE  LP 

C**  Generate  lineprinter  output 
INCLUDE  'BURGER. INC' 


*• 


DO  200  N=1 ,NMAX,NSKIP 
WRITE (LUO, 6000)  TITLE 
6000  FORMAT ( ' 1' / 

$  '  ' A60) 

VRITE( LUO, 6010)  T(N) ,  DT,  DX,  VIS,  EPSILON 
6010  FORMAT (T5' t  -'1PE11.3',  Dt  ='E11.3',  DX  ='E11.3 
$  ',  Mu  ='E11. 3' ,  Epsilon  ='E11.3) 

WRITE ( LUO, 6020) 

6020  F0RMAT(T5 ' I ' T9 ' X' T20 ' Ua ' T31 ' U' T42 ' %Er ror ' 

$  T53'ADR1'T64'ADR2'T75'ADR3' ) 

DO  100  I=1,IMAX,ISKIP 
PHIMIN  =  PHI SI 
PHIMAX  =  PHIS2 

PCTERR  =  100.  *  (  PHI(I,N)  -  PHIA(I,N)  )  /  (  PHIMAX  -  PHIMIN  ) 
WRITE (LUO, 6030)  I,  X(I),  PHIA(I,N),  PHI(I,N),  PCTERR 
$  ,  ADR1(I,N) ,  ADR2(I,N) ,  ADR3(I,N) 

6030  FORMAT ( 15 , 1P7E1 1 . 3 ) 

100  CONTINUE 

WRITE (LUO, 6040)  AMERR(N) ,  SIGMA(N) ,  RMSERR(N) 

$  ,  SUMERR(N) ,  SUMADRl(N) 

$  ,  PEAKERR(N),  PEAKADRl(N) 

6040  F0RMAT(T5' Arithmetic  mean  error  ='1PE11.3 

$  ',  Standard  deviation  ='E11.3',  RMS  error  ='E11.3/ 

$  T5' Integrated  ABS(error)  ='E11.3 

$  ',  Integrated  ADR1  ='E11.3/ 

S  T5'Peak  ABS(error)  ='E11.3 

$  ' ,  Peak  ADR1  ='E11.3) 

200  CONTINUE 

WRITE (6, 6050) 

6050  FORMAT ( ' 1 ' ) 

RETURN 

END 
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SUBROUTINE  MAC 

C**  Solve  using  the  MacCorraack  explicit  method 
INCLUDE  'BURGER. INC' 

DO  600  N=1,NMAX-1 

C**  First  check  stability  criteria 
DO  100  1=1, IMAX 

A(I)  =  C  +  B  *  PHI(I,N) 

CRITERIA  -  DX**2  /  (  ABS(A(I))  *  DX  +  2.  *  VIS  ) 

IF  (  DT  .GT.  CRITERIA  )  THEN 
C**  Vrite  a  warning  message  and  halt  execution 

WRITE (LUO, 6000)  DT,  CRITERIA,  N,  I,  PHI(I,N) 

6000  F0RMAT('  Dt  exceeds  stability  criteria'/ 

$  '  Dt  ='1PE11.3' ,  criteria  ='E11.3 

$  ' ,  N  =' 13' ,  I  =' 13' ,  U  ='E11.3) 

STOP  '  Program  execution  stopped' 

ENDIF 

100  CONTINUE 

C**  Predictor  step 

DO  200  1=1 , IMAX 

F(I)  =  C  *  PHI(I,N)  +  B  *  PHI(I,N)**2  /  2. 

200  CONTINUE 

C**  Set  boundary  condition  at  I  =  1  (Dirichlet  condition) 

PHIB( 1 ,N+1 )  .  PHI SI 
DO  300  1=2, IMAX- 1 
C**  Solve  interior  points 

FLUXC  =  F(I+1 )  -  F(I) 

FLUXV  =  VIS  *  (  PHI(I+1 ,N)  -  PHI(I,N)  )  /  DX 

$  -  VIS  *  (  PHI(I,N)  -  PHI(I-1,N)  )  /  DX 

FLUXS  =  VISA(I+1 ,N)  *  (  PHI(I+1 ,N)  -  PHI(I,N)  )  /  DX 

$  -  VISA(I,N)  *  (  PHI(I,N)  -  PHI(I-1,N)  )  /  DX 

PHIB(I,N+1)  =  PHI(I,N) 

$  -  DT  /  DX  *  (  FLUXC  -  FLUXV  -  FLUXS  ) 

ADR1(I ,N+1)  =  ABS  (  FLUXS  ) 

$  /  (  ABS(  FLUXC  -  FLUXV  )  +  ABS(  FLUXS  )  +  l.E-24  ) 

ADR2(I,N+1)  =  ABS(  FLUXS  ) 

$  /  (  ABS(  FLUXC  )  +  ABS(  FLUXV  )  +  ABS(  FLUXS  ) 

$  +  l.E-24  ) 

ADR3(I ,N+1)  =  ABS  (  FLUXS  ) 

$  /  (  ABS(  FLUXC  )  +  ABS(  FLUXS  )  +  l.E-24  ) 

300  CONTINUE 

C**  Set  boundary  condition  at  I  =  IMAX  (Dirchlet  condition) 

PHIB(IMAX,N+1)  =  PHIS2 

C**  Corrector  step 

DO  400  1=1, IMAX 

F(I)  =  C  *  PHIB(I ,N+1)  +  B  *  PHIB(I,N+1)**2  /  2. 

400  CONTINUE 

C**  Set  boundary  condition  at  I  =  1  (Dirichlet  condition) 

PHI(1,N+1)  =  PHIS1 

DO  500  1=2, IMAX- 1 
C**  Solve  interior  points 

FLUXC  =  F(I)  -  F(I-l) 

FLUXV  =  VIS  *  (  PHIB(I+1,N+1)  -  PHIB(I,N+1)  )  /  DX 


o  n 


$  -  VIS  *  (  PHIB(I ,N+1)  -  PHIB(I-1 ,N+1 )  )  /  DX 

FLUXS  «  VISB(I,N+1)  *  (  PHIB(I+1,N+1)  -  PHIB(I,N+1)  )  /  DX 

$  -  VISB(I-1 ,N+1)  *  (  PHIB(I,N+1)  -  PHIB(I-1,N+1)  )  /  DX 

PHI(I ,N+1)  =  0.5  *  (  PHI(I,N)  +  PHIB(I,N+1) 

$  -  DT  /  DX  *  (  FLUXC  -  FLUXV  -  FLUXS  )) 

ADR3(I,N+1)  =  ABS  (  FLUXS  ) 

$  /  (  ABS(  FLUXC  -  FLUXV  )  +  ABS(  FLUXS  )  +  l.E-24  ) 

500  CONTINUE 

C**  Set  boundary  condition  at  I  =  IMAX  (Dirchlet  condition) 
PHI(IMAX,N+1)  =  PHIS2 
600  CONTINUE 
RETURN 
END 
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SUBROUTINE  STATS 

C**  Compute  statistical  error  measures 
INCLUDE  'BURGER. INC' 

DIMENSION  ERROR (MAXI) ,  TMPl(MAXI),  TMP2(MAXI) 

DO  400  N= 1 , NMAX , NSKI P 

C**  Compute  arithmetic  mean  and  RMS  error 

SUM1  =  0. 

SUM2  -  0. 

DO  100  1=1 , IMAX 

ERROR(I)  =  PHI(I,N)  -  PHIA(I,N) 

SUM1  =  SUM1  +  ERROR(I) 

SUM2  =  SUM  2  +  ERR0R(I)**2 
100  CONTINUE 

AMERR(N)  =  SUM1  /  IMAX 
RMSERR(N)  «  SQRT(  SUM2  )  /  IMAX 
C**  Compute  variance  and  standard  deviation  of  the  error 

SUM3  =  0. 

DO  200  1=1, IMAX 

SUM 3  =  SUM 3  +  (  ERROR(I)  -  aMERR(N)  )**2 
200  CONTINUE 

VAR  =  SUM 3  /  IMAX 
SIGMA(N)  =  SQRT(  VAR  ) 

C**  Compute  the  integrals  of  the  absolute  error  and  ADR 
C**  Compute  the  peak  values  of  absolute  error  and  ADR 
PEAKERR(N)  =  0. 

PEAKADRl(N)  «  0. 

PEAKADR2(N)  =  0. 

PEAKADR3(N)  =  0. 

DO  300  1=1, IMAX 

ERROR(I)  =  ABS(  ERROR(I)  ) 

TMP1(I)  =  ABS(  PHI(I,N)  ) 

PEAKERR(N)  =  MAX(  ERROR(I) ,  PEAKERR(N)  ) 
PEAKADRl(N)  =  MAX(  ADR1(I,N),  PEAKADRl(N)  ) 
PEAKADR2(N)  =  MAX(  ADR2(I,N),  PEAKADR2(N)  ) 
PEAKADR3(N)  =  MAX(  ADR3(I,N),  PEAKADR3(N)  ) 

300  CONTINUE 

CALL  TRAPZ ( ERROR ( 1 ) ,  IMAX,  DX,  SUMERR(N) ) 

CALL  TRAPZ(TMP1( 1) ,  IMAX,  DX,  SUMPHI(N)) 

CALL  TRAPZ (ADR1(1,N) ,  IMAX,  DX,  SUMADRl(N) ) 

CALL  TRAPZ(ADR2(1,N),  IMAX,  DX,  SUMADR2(N)) 

CALL  TRAPZ (ADR3(1,N) ,  IMAX,  DX,  SUMADR3(N) ) 

400  CONTINUE 
RETURN 
END 
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SUBROUTINE  TRAPZ(F,  IMAX,  DX,  SUM) 

C**  Perform  trapezoidal  rule  integration  of  a  function  defined 
C**  by  a  table  of  equispaced  values. 

C**  Parameters: 

C**  F  Array  of  values  of  the  function 

O*  IMAX  Number  of  values 

C**  DX  Uniform  spacing  between  values 

C**  SUM  Estimate  of  the  integral 

IMPLICIT  REAL*8  <A-H,0-Z) 

DIMENSION  F(IMAX) 

SUM  =  0. 

DO  100  1=1, IMAX- 1 

SUM  =  SUM  +  DX  /  2.  *  (  F(I)  +  F(I+1)  ) 

100  CONTINUE 
RETURN 
END 
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SUBROUTINE  TRIDIAG(IL,IU,B,D,A,C) 

C**  Tri-diagonal  matrix  solver  (using  the  Thomas  algorithm) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  A(IU-IL+1) ,  B(IU-IL+1 ) ,  C(IU-IL+1),  D(IU-IL+1) 

C**  Reference: 

C**  Anderson,  DA,  JC  Tannehill,  and  RH  Fletcher.  "Computaional 
C**  Fluid  Mechanics  and  Heat  Transfer."  McGraw-Hill  Book  Company, 
C**  1984.  pp.  549-550. 

C**  IL  «  subscript  of  first  equation 

C**  IU  »  subscript  of  last  equation 

C**  B  *  coefficient  behind  (to  left  of)  diagonal 

C**  D  -  coefficient  on  diagonal 

C**  A  *  coefficient  ahead  (to  right  of)  diagonal 

C**  C  *  element  of  constant  vector  (right  hand  side) 

C**  Establish  upper  triangular  matrix 
LP  =  IL  +  1 
DO  100  I=LP,IU 
R  ==  B(I)  /  D(I-l) 

D(I)  =  D(I)  -  R  *  A(I-l) 

C(I)  =  C(I)  -  R  *  C(I-l) 

100  CONTINUE 

C**  Back  substitution 

C(IU)  =  C(IU)  /  D(IU) 

DO  200  I=LP ,IU 
J  =  IU  -  I  +  IL 

C(J)  a  (  C(J)  -  A(J)  *  C(J+1)  )  /  D(J) 

200  CONTINUE 

C**  Solution  has  now  been  stored  in  C 
RETURN 
END 
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FUNCTION  VISA(I,N) 

Compute  artificial  viscosity  due  to  fourth  order  smoothing 

INCLUDE  'BURGER. INC' 

IF(  I  .EQ.  1  )  THEN 

Use  first  order  forward  difference 
PIM1  =  PHI(I,N) 

PI  =  PHI(I+1,N) 

PIP1  =  PHI(I+2,N) 

ELSE  IF(  I  .GT.  1  .AND.  I  .LT.  IMAX  )  THEN 
Use  second  order  centered  difference 
PIM1  -  PHI(I-1,N) 

PI  -  PHI(I,N) 

PIP1  =  PHI(I+1,N) 

ELSE  IF(  I  .EQ.  IMAX  )  THEN 

Use  first  order  backward  difference 
PIM1  =  PHI(I-2,N) 

PI  =  PHI(I-1,N) 

PIP1  =  PHI(I,N) 

ENDIF 
R  =  1. 

R  »  (  ABS(  PHI(I,N)  )  +  C  )  /  (  PIP1  +  2.  *  PI  +  PIM1  ) 
VISA  -  EPSILON  *  DX  *  R  *  ABS<  PIP1  -  2.  *  PI  +  PIM1  ) 
RETURN 


FUNCTION  VISB(I,N) 

C**  Compute  artificial  viscosity  due  to  fourth  order  smoothing 
INCLUDE  'BURGER. INC' 

IF(  I  .EQ.  1  )  THEN 

C**  Use  first  order  forward  difference 
PIM1  -  PHIB(I ,N) 

PI  -  PHIB(I+1,N) 

PIP1  =  PHIB(I+2,N) 

ELSE  IF(  I  .GT.  1  .AND.  I  .LT.  IMAX  )  THEN 

C**  Use  second  order  centered  difference 
PIM1  =  PHIB(I-1,N) 

PI  =  PHIB(I,N) 

PIP1  =  PHIB(I+1 , N) 

ELSE  IF(  I  .EQ.  IMAX  )  THEN 

C**  Use  first  order  backward  difference 
PIM1  -  PHIB(I-2,N) 

PI  -  PHIB(I-1,N) 

PIP1  =  PHIB(I,N) 

ENDIF 
R  *  1. 

C  R  =  (  ABS(  PHIB(I,N)  )  +  C  )  /  (  PIP1  +2.  *  PI  +  PIM1  ) 

VI SB  =  EPSILON  *  DX  *  R  *  ABS(  PIP1  -  2.  *  PI  +  PIM1  ) 
RETURN 
END 
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